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been employed to solve the coupled stress balance, mass balance 
and energy balance equations. Since the governing equations for 
the visco plastic deformation of a material ar;e similar to those 
of a non-Newtonian, incompressible fluid, the flow approach has 
been used for solution purposes. The velocity, pressure and 
temperature fields thus predicted are further processed to obtain 
the isothermal patterns, deformation patterns and the cutting 
forces which are of practical importance. Also, a parametric 
study has been conducted to investigate, the effects of various 
process parameters such as depth of cut, cutting velocity and rake 
angle on the temperature and deformation patterns in the plastic 
zone. One of the significant contributions of this analysis is 
that important factors such as thermal softening of material, 
variation of material properties with temperature etc. have been 
considered. A complete software with pre-processing, FEM analysis 
and past -processing has been developed as a part of the study. A 
number of experiments for force and temperature measurements in 
orthogonal metal cutting have been conducted. Using the results 
of these experiments as well as some results of the earlier 
researchers, the theoretical model has been validated. Also, the 
present theoretical model confirms the energy minimization 
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principle proposed by Merchant for estimating the chip thickness 
or shear angle. 
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CHAPTER I 


INTRODUCTION 

! 

1.1 GENERAL BACKGROUND 

Metal cutting forms the basis of a wide variety of 
machining processesemployed in modern engineering industry. Nearly 
half of the engineering products today are produced through 
machining at one stage or the other. In this age of wide spread 
automation, we are witnessing the trend that the demand for 
productivity and process optimization is increasing day by day. 
In view of the crucial role played by machining among all the 
mater ial -processing operations, a scientific understanding of the 
metal cutting process is vital for modern engineering technology 
and practice. 

A detailed analysis of the mechanics of metal-cutting 
which sheds light on the underlying plastic deformation processes 
and thermal phenomena, is helpful in the prediction of some of the 
input parameters such as cutting forces. Evaluation of cutting 
forces is of paramount importance for the following reasons : 

(i) to estimate the power requirement of a machine tool 

(ii) to select appropriate bearings, jigs and fixtures which 

can withstand the tool forces 

(iii) to survey the characteristics of new work and tool 

materials by investigating tool forces. 

' I 

Further, minimization of the power input for a certain 
metal removal rate CMRR) often forms the objective of optimizing 
the metal cutting operation. Also, an insight into the effects of 
various process variables like the tool geometry, tool material 


2 


properties, cutt ing ' velocity and temperature distribution neai' 
tool -tip region is useful in understanding processes such as the 

I 

formation of built-up edge, tool wear and tool vibration. 
Minimizing such undesirable effects can increase productivity and 
lower machines costs. 

Over a past few decades, a large volume of experimental 
data has been collected on various aspects of metal cutting, but a 
comprehensive theoretical understanding of this highly complex 
process is yet to emerge. In order to estimate some of the 
important parameters such as cutting force, one has to rely on 
empirical relations based on experimental measurements and 
material properties. 

The present work, in its right earnest, is a small 
effort to fill some existing gaps in the theory of metal cutting. 
A coupled viscoplastic analysis of the deformation and thermal 
processes in the vicinity of the tool-tip has been attempted here, 
to provide a reasonably accurate estimate of the most important 
input parameters such as cutting forces and the process variables 
such as chip thickness ratio and tool-tip temperature. The 
theoretical predictions have been compared against experimental 
results obtained as a part of the present study and also with 
those of earlier researchers, to validate the theoretical model. 

1.2 REVIEW OF PREVIOUS WORK 

' I 

The earliest systematic study of metal cutting process, 
dates back to the end of 19th century, when pioneering researchers 
like Mallock, Tresca etc attempted to throw some light on the^ 
mechanics of machining. However, a major contribution to this 


area which has an iti/portant influence on the practical development 
of machining was possible only when Taylor CIU published his work 
in 1907. He was particularly concerned with tool wear and tool 
life and seems to have been one of the first to recognize the 
influence of temperature on tool wear. He developed a 

relationship between tool life and speed of machining which is 
still in use today. 

After this work, the major contribution in the area of 
metal cutting commenced from the post war period. The first 
attempts towards developing an understanding of the mechanics of 
metal cutting were performed by analysing the chip-formation 
process. As early as 1945, a pioneering work by. Merchant C21 and 
Piispanen C33 was presented on the so called classical "single 
shear plane model". This model assumes that the chip is formed 
due to plastic deformation occurring along a single plane, known 
as the shear plane. The most advantageous aspect of this model is 
that it enables the determination of the average yield shear 
stress and slip velocity along the shear plane, purely through 
simple geometric constructions. It also established beyond any 
doubt that metal cutting is basically a shear deformation process. 
The Merchant -Pi ispanen model was subsequently improved by Palmer & 
Oxley C43 and Okushima S; Hitomi C53. These studies suggested that 
though the formation of chip occurs due to flow of the metal under, 
shear, the deformation is not limited to a single plane; it occurs 
in a narrow region approximated by two parallel' planes. This 
model later came to be known as the thick shear zone model. In 
view of the small thickness of the shear zone in comparison to its 
length, this model assumes that the stake of stress within the 
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deformation region is uniform simple shear. 

Around the same period of the above mentioned 
theoretical models, efforts were going on to validate these models 
with experimentation. Keceloglu C6] was the first person to 
observe the process through photo-micrographs in 1958. The 
plastic zone where the chip formation takes place was photographed 
at short intervals using a quick stop device. By noting the 
deformation of the grain boundaries, the size and the shape of 
plastic deformation zone were estimated. This study established 
that the plastic region could be approximately represented by a 
thin par al le 1 -sided zone. Later in 1969, a more effective way of 
using photo-microgr aph technique was proposed by Stevenson and 
Oxley C7]. Assuming a thin par al le 1 -sided zone, they observed 
pr^inted grids on the material before and after deformation. From 

I 

the streamlines of metal flow, the strains and strainrates were 
calculated. The important process variable of flow stress Cor 
effective yield stress) could also be evaluated as a function 9 f 
strain by this method, for a wide range of strainrates. Later 
investigators [8,91 proposed that the flow stress is a strong 
function of strainrate as well. 

J.T. Black CSl proposed an approach based on the theory 
of dislocation mechanics. He suggested that the flow stress is 
composed of two parts i.e. thermal and non-thermal. The 
observations of the study indicated that the plastic zone is 
divided into alternate shear and lamella bands. The shear bands 
(or shear fronts) give rise to the dominant non-thermal part while 
the lamella regions lead to a thermal part of the flow stress- 
Evaluation of flow stress was based on the stacking fault energy 


o 


(^SFE). This approachr mostly used by metallurgists, is not very 

I 

useful to the continuum mechanics based theories. 

Von Turkovich C9, 103 first observed peculiar form of 
plastic flow, extremely localised in a small region in the ctiip 
ahead of the tool edge along the rake face. This was called the 
secondary shear deformation zone. The experimental studies 
conducted so far, have established that material behaves like a 
highly viscous fluid within the deformation zones and as a perfect 
iolid (with infinite viscocity and zero strain rate) outside them. 
The viscocity of such a fluid, depends strongly on material 
deformation rate. 

The credit for the effective usage of the FEM technique 
to viscoplastic analysis goes to Zienkiewicz ct.al. Cll, 12, 133. 
In 1973, these researchers applied FEM for the modeling of metal 
forming and extrusion problems. In metal cutting, the FEM 
technique was employed by Tay and Stevenson [14, 153 for 
predicting the temperature field for orthogonal cutting. They 
used hyperbolic streamlines in the primary deformation zone for 
calculating velocity fields and strainrates. They also derived 
semi empirical relations for computing heat generation. However, 
these studies do not provide a coupled analysis of plastic 
deformation and thermal processes during metal cutting. 

Muraka and Hinduja C163 predicted the temperature field 
variation during metal cutting and verified the average 
temperature rise by experimental measurements. These authors 
obtained the temperature distribution near the tool tip for a wide 
range of cutting conditions. Balaji C173 applied the FEM 


technique to evaluate the temperature field^ within the work 
material and the tool, for coated carbide tools. He also, has 
provided the experimental verifications of the results obtained. 
Mallikarjun Sarma C183 used the FEM-based viscoplastic analysis t'o 
predict velocity, pressure and temperature fields. But, his work 
did not consider some of the major aspects such as thermal 
softening and the variation in the material properties such as 
thermal conductivity, specific heat, yield stress etc. Still, it 
g^^ve a reasonably good estimate for the temperature field and the 
size of the primary deformation zone. In the most recent times, 
Komvopowrs etal. C31] have analyzed the deformation processes in 
metal cutting using elasto-plastic approach. The deformation 
pattern obtained here seem to reasonably good but, the thermal 
aspects as well as the experimental validation and estimation of 
important quantities such as forces is not attempted here also. 

Most of the FEM solutions available at present on metal 
cutting do not provide a complete theoretical analysis developed 
from the first principles. It is also observed that hardly any 
FEM work speaks of the parameters of practical importance such as 
cutting forces. Also very often gross simplifications have been 
invoked which are necessitated by the complexities of the process 
as well as the domain. These simplifications obviously affect the 
accuracy of the solution. In the present work, some of the 
limitations of the existing FEM solutions have been eliminated. 

1.3 OBJECTIVES AND SCOPE OF PRESENT STUDY 

In the present study, an attempt has been made to fill 
the gap in the analysis of metal -cut t ing by solving the coupled 



flow and energy equations. Stresses in the cutting zone have been 


modeled by considering visco-plastic material behaviour during 
deformation. Accordingly, the metal flow beyond the elastic limit 

\ I ' 

is treated similar to that of a liquid of very high viscocity. 
Viscocity itself is described as a strong, non-linear function of 

the strain-rate and yield stress of the work material. The 

/ 

material is considered to be purely strain-rate sensitive, with no 
effects of wor k -hardening. It must also be noted that material 
properties such as specific heat, thermal conductivity and 
uniaxial yield stress vary with change of temperature. Metal 
cutting is bound to involve large temperature variations. Hence, 

I 

to have a realistic picture of the process, the variations in 
material properties with respect to temperature have been 
considered. In particular, the effect of thermal softening has 
been incorporated by considering the reduction in the value of 
yield stress with increase in temperature. 

It is well established from previous experiments that 
the deformation is restricted to a very small region around the 
tool tip. Therefore, a solution domain extending up to a few 
millimeters from tool-tip into the workpiece and the chip has been 
considered. Beyond this domain, the material is assumed to be 
perfectly rigid. The tool is also taken to be perfectly rigid and 
sharp with no built up edge. The tool cutting edge is assumed to 
be orthogonal to the cutting direction. Further, for generating 
the domain and finite element mesh easily, the geometry of the 
problem has been simplified in the following manner. The curved 

ft 

top surface of the chip is given by a hyperbola with asymptotic 
parallel to the feed direction and the rake face. The bottom 


surface on the Qt^lGr hand, is in perfect contact with the rake 
face upto the contact point and it becomes parallel to the top 
surface beyond the point of contact. The chip thickness is 
determined as a part of the analysis by invoking the minimum total 
energy pr inciple. The inlet and exit boundaries are taken to be 
orthogonal to the feed direction or the rake face. Though the 
ch'ip curling process makes the real boundaries df the problem very 
complex, the simplifications employed in the present study are 
appropriate, as the focus of this investigation is to perform 

I 

coupled stress and heat transfer analysis for a given domain. 

Anotfier important assumption invoked in the present 
study is that all the powers expended in plastic deformation is 
converted into heat. Though some residual stresses/strains are 
retained both in the workpiece and the chip, this assumption 

I 

simplifies computation to a great extent without sacrificing the 
accuracy of predictions, when compared with experimental results. 

The governing equations of flow and energy after the 
incorporation of the above mentioned assumptions, have been solved 
using the Finite Element Technique. The velocities, temperature 
and pressure distributions are obtained and utilized to identify 
the shapes of isotherms and for the prediction of cutting forces. 
Also from the strainrates obtained in the domain, a reasonably 
good idea of the dimensions of primary and secondary deformation 
zone is obtained. 

A special software for FEM grid generation has been 
developed in the present study, since grid generation for this 
complex geometry encounter«| in the metal cutting problem is quite 


a difficult task. Some of the advanced concepts such as the 
hyperbolic streamline model proposed by Tay and- Stevenson C15] 
trans-finite interpolation etc. have been used in the grid 
generation algorithm. To have a good overview of the predicted 
results, the analysis has been carried out for various cutting 
conditions viz. different depths of cut, cutting velocities and 
rake angles- A number of experiments on orthogonal cutting have 
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been performed to verify the results obtained from the theoretical 
analysis. 

After developing a complete package consisting of the 
finite element analysis and grid generation, the effect of 
variation in chip thickness on the total energy consumed during 
metal cutting has been studied. This can be viewed as an 
excellent check for the present theoretical approach, since 
Merchants Cl] theory uses the energy minimization principle while 
evaluating the chip thickness or shear angle- 


CHAPTER ^ 


THEORE7 1 CAL FORMUI . AT 1 ON 

I 

Metal cutting involves very complex deformation 
phenomena. Extremely localised asymmetric deformation occurs at 
exceedingly high strains and strainrates, in a minute domain of 
few millimeters size. The metal cutting action itself is a 
consequence of high compression followed by shearing of the metal. 
Consequently, any theoretical treatment of the metal cutting 
process is bound to possess certain relaxations and 
simplificat ions. 

A two dimensional orthogonal machining under dry cutting 
conditions has been considered for analysis. In fact, orthogonal 
cutting is a particular case of oblique cutting, with cutting edge 
of the tool perpendicular to the cutting velocity vector. The 
tool is assumed to be perfectly sharp, with no builtup edge. The 
work material is taken to be ductile, resulting in the formation 
of continuous chips. The tool and the portion of work material 
outside the solution domain on the other hand, are taken to be 
perfectly rigid without any deformation. 

2.1 Plastic doformation and chip formation in orthogonal 

machining 

In orthogonal metal cutting, the fact that the width of 
the chip is large compared to the its thickness, renders the 
problem two-dimensipnal . Although the deformation in metal cutting 

I 

is similar to those of the metal forming processes a 
distinguishing feature is that the deformed material, is separated 
from the parent workpiece in the form of a chip. Hence, it is of 


when the 


prime importance to analyze the chip formation process, 

theoretical analysis of metal cutting is done. 

1 

The uncut material in motion is first interrupted by the 
cutting tool, causing severe compressive stresses in the material. 
After sufficient compression, when the resulting shear stresses in 
the material attain the yield value, plastic flow of the material 
occurs in the direction of shear forces; leading to the formation 
of the chip. The outward or shearing movement of each successive 
deforming element is arrested by workhardening and the movement 
is transferred to the next element. 

It is well accepted now that the predominant mechanism 
of metal cutting is the shearing process occurring in the zone 
called "Shear Zone". A detailed study of metal cutting reveals 
that two such zones exist in the domain. 

i) Primary shear deformation zone (PSDZ) 
ii) Secondary shear deformation zone (SSDZ) 

These two zones have been shown in Fig. 12.1). The 
primary shear zone has been observed to be wedge shaped and its 
thickness mainly depends upon cutting velocity. It has been 

reported that an increase in cutting velocity causes the thickness 
of PSDZ to reduce. The thickness or mean width of PSDZ is an 
important quantity which gives an idea of the rate of deformation 
and the amount of heat generated. The narrower the region, the 
greater the rate of deformation. 

’ I 

The secondary shear deformation zone lies adjacent to 
the contact patch between the chip and the tool. Additional 
deformation due to the friction takes place fiere. Eacti layer io' 
SSDZ moves at a different speed because of shearing in a direction 
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parallel to the rake face. Although the region is small as 
compared with PSDZ, the amount of heat produced is believed to be 
of a 'similar magnitude. An interesting point to note is that, 
the highest temperature occurs on the interface between SSDZ and 
rake face. The temperature in SSDZ is higher than in PSDZ mainly 
because, the metal flow velocity is very low in SSDZ and also the 
material which enters SSDZ has already been heated through PSDZ. 
Some of the researchers have attempted to estimate the thicknesses 
of these zones theoretically and experimentally C27j} . 

2. jS Viscoplasticity Concepts 

I 

The study of the behaviour of a ductile material under 
tensile load reveals that beyond the elastic limit when the 
material yields plastically, the stress (some times called as flop 
stress) becomes a function of strain-rate. Some of the properties 
of the solid undergoing plastic deformation are akin to that of a 
Non-Newtonian Fluid in the sense that, the relationship between 

I i 

the shear stress and the shear slrainrate is non-linear. 

I The material incompressibility is often assumed in the 

analysis of plastic deformation processes and is quite reasonable 
since, yield occurs primarily due to shear. Therefore, the sum 
of all the diagonal components of str ain -tensor (which represents 
volumetric strain) is zero. 

i.e. e..=e..+ e„„ + = 0 (2.1) 

11 11 22 33 

Deviator ic stress and deviatoric strain or strainrate 
are of vital importance in theory of plasticity. These are 


defined as 
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( 2 . 2 ) 

(2.3) 


where 


e . . 

1 X 


volumetric strain 

e . . 

1 J 


strain tensor 

. 

1 J 

= 

deviatoric strain tensor 

<y . . 

1 J 


stress tensor 

^ii/3 


isotropic stress 

1 J 


deviatoric stress. 


Various t^^eories of plasticity developed so far, attempt 
to relate the deviatoric part of strain (or strainrate) with 
deviatoric stress. 

An extensively used yield criterion in the analysis of 

plastic deformation processes is due to Von liises. According to 

this criterion, the second invariant of deviatoric stress tensor 

(represented by J„) determines the criterion of yielding. 

2 

Mathematically, yield occurs when 



where, K is a material yield constant. 

In the plastic deformation, the total strain deviation 
consists of two parts i.e. elastic and plastic. But, in the 
problems such as metal cutting or metal forming involving huge 
plastic deformations, the elastic strain deviation makes a 
co'ntr ibut ion of just 2 to 3'/- only and can therefore be neglected- 
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Hence , 


e . 
1 J 


1 j 


C2.5) 


The total plastic strain e occur r ing 

1 J 

period of time t is given by the area under 
strain-rate te. .) and time tt>. Thus, 

i- J 


e'';’ = ( 0 ) + r 6. . di 

1 J 1 J J IJ 


for a loading 
the curve of 


( 2 . 6 ) 


C P ^ 

where e. . (0) is the plastic strain at the initial time t = 0. 

From the principle of material incompressibility (from eqs. 
2 . 1 - 2 . 5 ) 


e 


(p) 

i i 


= 0 


(2.7) 


' The strainrate tensor ef*^^ is thus, , same as 

ij ' ' ij 

Sine, only the rate of plastic strain is specified by above 


equations for a varying loading sequence, the total plastic strain 
C p) ) 

e; . can often be replaced by the algebraic sum of successive' 
1 J 

strain increments. 


The addition of all the mirmte strain increments to 
obtain the total strain forms the essence of the the incremental 
theory of plasticity. 


a. 3 Mathomatical Formulation 

2.3.1 Constitutive rolation for viscoplastic flow 

The metal undergoing plastic deformation can be 
described as a Von-Mises type of strain-rate sensitive, non-strain 
hardening, visco-plastic material. This, in turn, is similar to 
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an ^ incompressible, non -newton i an fluid for which the above 
mentioned viscoplasticity concepts can be easily applied. 

An important aspect to be noted here is that, the 
viscocity of such a "Fluid" is dependent on the current locals 
strainrates, total accumulated plastic strain undergone by the 
material and the local temperature. It has been assumed in the 
present analysis that work hardening is absent. Thus, for a 

' j 

purely strainrate sensitive material, the constitutive relation is 


or . = 2/Li er . (2.8) 

1 J X j 

where, fj = viscocity. 

er . = deviatoric part of strain rate tensor 
X J 

The strainrate tensor, for a two dimensional case can be written 
as 


^11 ~ 

(2.9a) 

^22 " dv/dy 

(2.9b) 

. - . fdu . dv T 

12 21 l^^y dx J 

(2.9c) 


where, u and v are the velocity components in x and y directions. 

The viscocity (J for a strainrate sensitive material can 
be written in the form of the constitutive relation given by 



■/3 e 


( 2 . 10 ) 


where a = yield stress of the material 
y 

r Z, T\ = Physical constants which define the viscoplastic 
characteristics of the material. 
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Effective yield stress which is often called as Flow 
stress (.O') is a function of the uniaxial yield stress and 
strain-rate Ce). Thus, 



From eqn. (2.10) and (2.11) we get. 


( 2 . 11 ) 


a - Vs fJB . ( 2 . 12 ) 

The strain rate invariant e can be expanded in terms of strain 
rate tensor as. 


e 



( 2 . 12 ) 


For a two dimensional case, above eqn., (2.12) reduces to. 


e 12 ^^^ 12 ^ 21 ^ 

2.3.2 (governing Equations : 

In the present analysis, the metal cutting problem is 
being viewed with the perspective of the viscoplastic deformation 
process. Especially, the velocity fields, pressure field and the 
teq-iperature field in the plastic domain need to be predicted. 

From the consideration of material balance, stress 
balance and heat balance, the steady state governing equations for 
the velocity, pressure and temperature fields are given by : 


(2. 13) 


(2. 14) 
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V. V = O (From mass balance) 


C2. 15) 


and 


pi ", V V) = V. <?■ + f (From stress balance) 


(2. 16) 


2 .(W.VT)+ Q = pQ (V .V T) (From energy balance) (2.17) 


where 


P 

y 

f 

k 

T 

Q 

C 


Density of the workpiece material 
The velocity vector 
Stress tensor 
Body Force tensor 

Thermal conductivity of wor k -mater ial 
Temperature 

Heat generation rate per unit volume. 
Specific heat of workpiece material. 


Flow Equations t 

The stress tensor in the momentum equation can be 
expressed as. 


n = - p;l + d 


(2. 17) 


where, p is the pressure and d is the deviator ic part of the 
stress tensor. 

The deviator ic stress tensor d can be written as. 


d = Htv y + y y') - ± (y. y) i 


(2. 18) 


where, V is the transpose of V V. 

A/ #w 

Simipl if y ing eqn. (2. 18) r using inconr^pressib i 1 i t y condition^ we 
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can write. 


§ 


T 

-PI +■ u ('7 V +- V V < ) 

^ r%» 


Thus, the momenturri equation gets the form. 


C2. 19) 


p (V.V V) = 9. |-P^ + fj (9 y + V y^) | t2.20> 

where, the body force vector f (which usually is the weight of the 
material per unit volume) has been absorbed in to the pressure. 

Component. Form of Flow equations 

Material balance equation for a two dimensional case can 
be written as. 


=0 

dx dy 

where, u and v are the velocity components in x and 

The stress tensor ^ dimensional 

expressed in component from as : 


( 2 . 21 ) 

y directions, 
case, can be 


a 


a n + O' 1 j + O' j 1 + O' jj 
XX xy yx yy 


(2.22) 


Thus, the component forms of the momentum equation can be written 
as follows : 


X-momentum 


and 


r ^u <?ui r ^ ^ ^ 

I U-s— I = < .=—(0' ) +■ (o’ ) >- 

dx 9y} Sy. XX dy yx J 


(2.23) 


Y-momentum 


r av ^ ^v’) r ^ , 

Ox OyJ I Ox xy 


Oy i y 


} 


(2.24) 



Stress tensor conjponents are given by : 


and 


O' 


v.^ 


p 



a = O' 
xy yx 


d\i 

dy 



o 

yy 


p 




£v 

dy 


(2.25) 


Using (2.25), the final forms of mome,ntum equations are 


P 


P 



(2.27) 

! 


(2.28) 


Energy Equation ' 

The objective of solving the energy equation is to 

! 

predict the temperature field. It is thus obvious that, the heat 
generation must be evaluated carefully- In metal cutting, it can 
be reasonably assumed that the heat generation is due to (i) 
plastic work in PSDZ, (ii) Friction in SSDZ. It has been assumed 

Co-\so Sow^flE. p\ a*Wc Woir 

that all of the work done is converted into hea,t , which seems to 
be quite reasonable because, the material density does not change 
considerably and also, the work hardening is not severe due to 
high temperature conditions. 


Volumetric heat generation rate can mathematically be 

r 


expressed as 
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where, t = Effective yield sheo-r stress, 
and e = strain rate invariant. 

Using Von-liisses criterion, f.2.29> can be simplified as, 

Q = — - e (2.30) 

yi 

Thus, using (2.30) and (2.17), the final form of energy 
equation becomes. 


p^T ^ a T ■ 

. O' “ ^ 

r aT ^ aT'i 

Ux^ av^ - 

+ — e = p C 

['^ax '^ayj 


(2.31) 


2.3.3 Non-Dimenfeionallzation t 

As long as the number of parameters involved in a 
problem are few in number, it is realistic to stay with the 
dimension al form of variables, which enables ready comparison of 
predicted results with real values. However, when the number of 
variables is large, analysing the influence of each parameter on 
the final result becomes a tedious task. The subject of metal 
cutting is known to involve parameters, which exhibit very complek 
interdependences among themselves. To circumvent such hardships 
and facilitate a detailed parametric study, it is desirable to 
non-dimensional ize the variables of the problem with suitable 

scaling factors. 

'' 

A real difficulty which often arises during 
non-dimensional ization, is the proper choice of a scaling factor 
for particular quantity. Though some general guidelines do exist^. 


it can be stated that. 

the 

primary guide 

in this 

pr OCHBS 

is 

intuition alone, since 

there is no 

unique 

way 

for 

non -dimensional izat ion. 

However r the best 

possible 

way is 

to 



critical examination and 


give, the physics of the problem, a 
incorporate as much of tlie physical char ac ter ist ics of the problem 
as possible while selecting the scaling factors. 

For the present problem, to begin with, the velocity 

components u and v are non -d imensi onal ized with the cutting 

velocity (V); the distances with depth of cut (tl? and the 

pressure with the uniaxial yield stress at ambient temperature 

(c ). The scaling factors for the remaining variables will be 
yo 

selected conveniently, during the course of nori-dimensional izat ion 
itself. 

Thus, ws define 

= x/i; = y/t? = u/V? = v/V 

where ^ denotes a dimensionless quantity. Using the above 
definition in material balance equation^ 



Momontum Equations 
X -Momentum equation 

x-moment urn equation in dimensional form is given by 



wherer P = Density of the materil 
jj rr viscocity 

Atr this stager the viscocity can be scaled as shown 

1 

below : 




i<y t/V) 

YO 


o 4- 

y 


f ^ — 1 

^ -/3 r J 


1 /n 


-/3 


ia t/V) 
yo 


or 


where. 


(o- /<y ) + -J I (e*)^ 

^ ^yo 1 ^ i 

V 3 e"" 


/n 


» 


O' ^ ^ 1 /n 

y + B (e ) 


ya e* 


(2.34) 


a 


a /a 
y yo 


B 


= _ 1 _ [A JL_\ 

^yo 1 y3 


1/n 


constant for the prescribed 
material and cutting conditions 


and e*= e/(V/t) 

Now, non-dimensional izing each quantity of the equation t2.33) we 
get 


pV 


2 f »du . »3u 




+ V 


dx 


dy 


f 5 4.^ r 

I dx^ dx^ dy^ L 




» [5u . erv 


dy” dx 


:]]} 


C2-35) 


or^ r 


/ 





I 

a 






d>{ 





a ^ 

ax'" 



(2.36) 

2 

where, a - a / (pV ). 
yo 


Y-Momentum equation 

Similarly, on the above lines the y-momentum equation 
can be reduced as. 




d>i 




•¥ V 




J = f;« 


(-p«+2p* ^^) 
9y^ dy^ 


(2-37) 


Energy equation ^ 

» 

The dimensional form of energy equation in the present 

study is. 


I 

where. 


a 

dy. 


(k 


dl 

di< 


) + 


a 

ay 


(k 


aj 

ay 


) +Q=pC (u 


aT 

ay 


+ V 


fl' 

ay' 


(2.38) 


k = Thermal conductivity of the work-material 

C = Specific heat of the work material 
P 

0 = Heat generated. 


The 

be done using 


normal izat ion of temperature T in this equation 

some reference temperature T _, as 

y s T 


can 
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.% T 


ref 


can be conveniently chosen as shown below. 
Non-dimensionalising each quantity of euqation (.2.38), we get 


r 

■^ref { 


— ( k — ) + — ( k — ) + T- O’ — e 

dx* ax* ay* ay* ^ VS 


} 


pVC T 
po re 




» ai* ^ w ai* 

ax* ay 


:>] 


(2-39) 


Multiplying above equation by t /(pC }y 

po 


^ 't ( 

3C ref 1 

DO V 


. ^ a_ ar , U vt , I ?_ 

^ dx** ax” ay” ay” J ^po 


= Vt T 


re 


f ['^P 


^ ^ • 
* ^ V* II ) 

ax* ay* ■ 


(2.40) 


Taking, <y /ip C ) = Tref 
yo po 

and k / (p C ) = a (Thermal dif f usivity ) , 
o po 


the final form of dimensionless heat balance euqation 


becpmes 


{ 


a aT*. ^ a ,,s ^t*,) ^ d 

ax* ax* ay* ay* / ^ Vs 




» aT* ^ jf aT* 

4. V 

0x dy 


,] .0 


(2,4i3 


wherGy 


Non --d imensional thermal conductivity - k/k 


Non-dimensional specific heat = C /C 

p po 


Peclet Number = Vt/a 




Fig. 2-3 Concept of hyperjjollc streamline used in generating 
VIII surface of the solution domain 
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VII. 

The free surface of the chip, where the complex chip 
curling effect with serrations is observed, liow been modelled using 
the Hyperbolic streamline concept used by Tay et al. tl53. This 

forms the surface VIII. 

Hyperbolic stream line approach 

This was used by Tay, Stevenson et al. (15) to find the 
distribution of shear strain rates in the deformation zones. The 
eighth surface of the present domain, very closely approximates a 
hyperbolic streamline used in the above work. 

The hyperbolic streamline is of the form 

tan r — xy = a (2.6.1) 

which is illustrated in Fig. (2.3). With this equation, the 
origi'n of the Cartesian Co-ordinate system which lies on AB is 
located at the intersection of the two asymptotes one of which is 
parallel to the cutting direction while other is inclined at an 
angle (90° - y) to the cutting direction. 

Thus, the equation for the hyperbola defined in Fig. 
(2.3) is given by the equation, 

(y-mx-y) (y-1.0) = (yl) - mxl - y) (yl - 1.0) (2.6.2) 

and slope 


dy 


= m(y-1.0)/C(y-1.0) +(y-mx-y)3 
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k ?: C are the thermal conductivity and specific heat of the 
o po 

work material at ambient temperature respectively. 

I 

2. 4 Selection of Solution Domain for Viscoplastic Analysis 

It is well established from previous experiments that 
the deformation is restricted to a very small region around the 
tool tip. Therefore, a solution domain extending up to a few 
millimeters from thetool tip into the workpiece and the chip has 
been considered. Beyond this domain, the material is assumed to 
be perfectly rigid. The tool is also taken to be perfectly rigid 
and sharp with no builtup edge. The tool edge is assumed to be 
orthogonal to the cutting direction. Further, for generating the 
domain and the finite element mesh easily, the boundary of the 
domain has been divided in to eight surfaces, as shown in Fig. 
(2.2). Surface I is the entry section where material enters the 
plastic solution domain and is considered to be plane, 
perpendicular to cutting velocity. Surface II, is selected in a 
plane parallel to the cutting velocity and restricts the boundary 
of the plastic zone below depth of cut. Material leaves the 
plastic zone at surface III which is parallel to the surface I. 
Surface IV is formed by the machined surface and is exposed to the 
atmosphere due to the the clearance angle of the tool. 

The rake surface, along which perfect sticking contact 
exists between the chip and the tool form the V surface 
(sticking contact region) and rest of the contact region i.e. 
slipping contact region forms the surface VI. 

The chip is considered to be bounded by plane parallel 
surfaces beyond contact point? which is represented by Surface 
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where, m = 


s i n ( 90° -r ) 
cos (90° -7-) 


y = L y - ay 1 - ( >: -a>; 1 ) m 3 / ( 1 -a ) 

o o 


(2.6.3) 


and a, x.ry., x , y etc. have been defined in the Fig. (2.3) 

1 1 o o 

1*! Estimation of Contact Length t 

In continuous chip formation process, the chip remains 
in contact with the rake face of the tool over a certain length. 
This contact length is known as natural chip-tool contact length. 
The importance of the chip tool contact length in the mechanics of 
continuous chip formation was not realized for a long time, 
because of the fact that most of the theories of metal cutting, 
such as those developed by Merchant, Lee and Shaffer etc. did not 
include this as a parameter. But it has since been established 
that, the force of friction acting on the chip-tool interface 
depends upon the length of chip-tool contact. 

/ 

Some amount of research is available on the contact 
length. Venkatesh et al. (20) have measured the contact length 
successfully using the Wheatstone bridge method. Almost all th^ 
theoretical analyses for the determination of contact length in 
natural contact cutting are based on the equilibrium of chip under 
the action of moments of forces acting on the tool -rate face and 

I , i 

on the shear plane. Hahn (29) has developed a relation for the 
determination of the natural chip-tool contact length, making use 
of merchants theory as well. The values predicted by this 
relation : have been found to matcfi with experimental values 
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reasonably, especially for the small depth of cuts of the order of 
-1 -2 . 

10 to 10 millimeters. The relation is given by, 

, _ t sin 

' c sTh'^ cos (3 

= Natural chip-tool contact length i 

= Depth of cut 

= Shear angle 

= Friction angle 

= Rake angle 

In the present analysis, chip-tool contact length is an input 

parameter. 

A simple method such as applying lamp black on the rake 

I 

surface and then measuring the scar produced by the chip on it 

under travelling microscope, may be used for measurement of 

contact length- Some such measurements have been made in the 
present study also. 

In the solution domain, surface V eVctends over the 

sticking contact length and surface VI over the slipping part of 

contact length. In the present work, both sticking and slipping 
contact length are assumed to be equal. 

2.5 Boundary Conditions 

The non-dimensikonal form of the boundary conditions are 
discussed below. 

Refer ing to Fig. 2.2, for convective heat loss on boundaries IV 
and VIII, the dimensional for of boundary condi-feion is 


whet^e, 

t 

<P 

0 

a. 
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~k 

= h 

dn 

( T-T ^) 

(2.42) 

where h = Heat transfer coefficient. 


Non dimensional i zing 

1 

each quantity of eqn (2.42), 

1 


-k 

« £T* 
an * 

h« ( t^-T^ . ) 
amb 

(2.43) 

where 





= 

T/T - 
ref 


h* 


h t/k 

o 


T^ = 

amb 

^amb ^ ref 



On ithe chip tool interface (surface V and VI), the frictional heat 
generated (Qp) is assumed to be partitioned according to the 
thlirinal conductivity values of the chip and the tool. The leads 
to the boundary conditions of the farm 


-k 


£I 

dri 


(Wt] 




(2.44) 


where Op = Frictional power dissipation 

= Tp X 

k, k^ = thermal conductivities of the chip and the tool 
materials 

T_ = Effective frictional shear stress (= — or mpe) 


Tangential velocity of slip 


— Frictional heat partition ratio 



1 
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Non -d imensional i 2 ing each quantity of (eqn. (2.44) 


, , » ref dl w 

■k k — T — - <p. y. a x t_ x V x V,. . 

o t_»^kyoF T 

an 


Simplification of above equation gives, 


-k — = <p P T_ V- 

- » k e F T 
dn 


(2.45) 


where p = peclet number |= Vt/(k /pC ) 
e L o po J 


On surfaces I, II and III 


lere, . = T . /T 

amb amb' ref 


At the far end of the chip (surface VII), 


(2.46) 


(2.47) 


This condition implies that, the temperature variation in the chip 
flow direction is negligible and to the absence of heat generation 
or any major heat loss mechanism. 

. The dimnsionless velocity conditions are ; 

' u* = 1, v** = 0 on srufaces I, II, and III 

V* = 0 for surfaces IV, V, VI and VIII 
n 

The chip velocity V^ also can be non-dimensionalized as 


resulting in 
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u 

= 

c 

cos 

0 ( 

V 

If 

<: 

n 

sin 

o< 


For surface VII 


(2.48) 


Frictional stress condition on Surface V becomes. 

1 1 /n 


F O' 


yo 
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-/3 e 
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1/n 
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f * 
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+ (e*) 


l/n| 


(2.49) 


where B* 



I 

(2.50) 


O' 

Y 

m 


a /a 


y yo 

Friction factor. 


^ The dimensionless governing equation and boundary 
conditions formulated above have been solved by the application of 
the Finite Element Method. 


a. 6 


Energy Minimization Principle 


Chip thickness is a significant parameter of 
solution domain. With the change in chip thickness or 
thickness ratio, the size of the solution domain changes. 


the 

chip' 


According to Merchant (1), whose theory is based upon 
the "minimum energy principle^, the shear plane is located where 
the least energy is required for shear- This principle of 
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Fig. 2-6 Variation of Uniax Yield stress with temperature 
(for plain carbon steel ) 
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minimum energy has been verified in the present analysisr by 
varying the chip thickness and finding the total plastic work and 
frictional work involved for each case. 

2. 7 Estlmat-ion of Material Properties 

I 

A careful study of the governing equation reveal that 
the material properties such as the uniaxial yield stress, thermal 
conductivity and specific heat affect the solution to a 
considerable extent. All the above mentioned properties vary with 
the temperature. In particular, at higher _ temperature the 
variations become quite significant. 

In the present analysis, the solution domain is a minute 
region surrounding tool tip totally under plastic deformation. As 
already mentioned elsewhere, the temperatures in this region are 
quite high. Thus, to have a realistic picture of the results, the 
variation in these properties at higher temperatures is 
incorporated in the. analysis. 

The variation in thermal conductivity and specific heat 
of mild steel with respect to temperature is shown in Fig (2.4) 
and (2.7). This is taken fftpm published data (30). The variation 
is uniaxial yield stress of mild steel with temperature is shown 
in Fig. (2.6) and has been taken from reference (28). 

Metal cutting being a process involving high strain 

3 

rates of the order of 10 /S, the variation in yield stress is 
selected accordingly. 

These variations are curve fitted and property 

correlations have been derived. For thermal conductivity and 

\ 

specific heat, piecewise equations of a straight tine of the type 
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y = mx + c have been fitted for different teiriper at ure ranges. For 
uniaxial yield stress, on the other hand a parabolic equation of 
the type 

y = 9.6991 — 2137677X + 1521.180 
has been derived using the least square method;, where 
X = temperature in 

2 

and y = uniaxial yield stress in N/m . 



CHAPTER 3 


I 


I 


FINITE ELEMENT ANALYSIS 


3.1 Introduction to Finite Element Technique 

The finite element technique, whose engineering birth 
and boom in the 1960s was due to the application of digital 
computers to structural analysis, has spread to a variety of 
engineering and physical science disciplines in the last decade. 
In many engineering problems, it is always not possible to obtain 
a closed form exact solution. Thus, for such difficult problems, 
it is possible to get the best approximate solution (i.e. very 
close to exact one), using the powerful numerical techniques such 
as finite element method. 

The basic concept of the finite element method is one of 
discretization. The finite element model is constructed in the 
following manner. A number of finite points are identified in the 
domain of the function, and the values of the function and its 
appropriate derivatives these points are selected as unknown 
variables. The points are called as nodes. The domain of the 
function is represented approximately by a finite collection of 
subdomains called finite elements. The domain is then an 
assemblage of elements connected together appropriately on their 
boundaries. The function is approximated locally within each 
element by continuous functions that are uniquely described in 
terms of nodal point values associated with the particular 
element. The path to the solution of a finite element problem 
consists of five specific steps : 



' (i) identification of the problem i 

( i i .> discretization and definition of the element 
( i i i lestabl ishment of the nodal equations for one element 
(iv) assemblage of elemental contributions to nodal equations 
and incorporating boundary conditions and 
Cv) the numerical solution of the global equations. 

The formation of i element equations is done in four ways. 
(1) direct approach. (2) variational method t3) method of 
weighted residuals and (4> energy balance approach. 

In the present analysis, GalArkin^'s weighted residual 
method has been used. 

3., 2 APPLICATION OF FEM 

3.2.1 Grid Generation : 

The first step in the finite element method is the 
selection of solution domain and its discretization. The choice 
of the solution domain boundaries for the metal cutting problem 
was discussed in previous chapter (see. Fig. 2.2). Grid 

generation for such a domain is quite a tedious task. 

Moreover, it can be seen that the variation in depth of 
cut, chip thickness, rake angle or contact length leads to a 
totally new mesh. Thus, it is obvious that, unless one has an 
efficient way of grid generation, the analysis can not be done 
conveniently. Thus, grid generation forms one of the important 
aspects of the present analysis. 

Considering the unsymmetry of the domain, it is divided 
into two zones, namely zone 1 and zone 2. The shapes of these 
zones are defined as shown in the Fig. (3.1). 
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of = Rake angle 
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Fig. 3-1 Division and size of the solution domain set for 
grid generation 
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Figi 3-2 Concept of trans-finite interpolations 
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Coordinate data and nodal connectivity of elements for 
each zone are calculated separately first and global elemental and 


nodal 

numbers for both zones put 

together are 

assigned 

later. 

To 

begin 

with. 

the coordinates for all the 

boundar y 

nodes 

are 

pr escr 

ibed 

suitably and then 

using a 

special 

type 

of 


interpolation, known as trans-finite interpolation technique, 
coordinate data of interior points have been generated. This 
technique is based on transforming an arbitrary 4-sided domain 
into a square domain in terms of body-fitting coordinates t? , y?) . 

The Cartesian and body fitting coordinates are shown in Fig. 3.2. 

The procedure for generation x and y values for a 
typical point H is illustrated in Fig. 3.2. The domain ABCD is 
mapped into a unit square Clxl) in the transformed domain. 

If H is the intersection point of the grid lines PQ and RS (which 
are ry = const. and ? = const. lines respectively), by 

bi-directional interpolation the x and y coordinates of H can be 
obtained in terms of the x and y coordinates of the boundary 

points P,0,R and S. For instance, if ? and 1 -? are the distances 

1 

of H from P and Q respectively in the transformed domain, then for 
one -d ir ect ional interpolation 

= ( 1 -?) Xp + ? Xg (3.1a)' 

= (1-?) y + ? X (3.1b) 

However, for two directional interpolations the interpolations 
must be applied in both ? and 77 directions suitably and 
ard then obtained in terms of the coordinates of the boundary 

points P,Q,R and S. Similarly for obtaining the coordinates of 
all the interior points, the appropriate boundary points (Fig. 
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Fig. 3-3 A sample of the mesh generated by the grid generation soft 
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3.1) 35 4--5idedlT by transfinits interpolation the interior 

coordinates can be generated. A typical mesh thus generated 
is shown in Fig. (3.3). 

In the present analysis, a grid with ^36 nodes and 111 
elements has been used. It was observed that, lesser number of 

elements and nodes and large sized elements, often led to the 

/ 

numerical difficulties such as negative determinant value etc. 
Provision has been kept in the program, to change the number of 
elements and the size of the elements to overcome such problems. 

; I 

3.2.2. DERIVATION OF FINITE ELEMENT EQUATIONS 

' The finite element equations for plastic flow of metal, 

are better derived from the basic vector equations so that the 
handling of the traction boundary conditions at the chip-tool 
interface could be easily implemented. 

FLOW EQUATIONS 

The vectors momentum equation is given by 

p(V.W) - 7. e - f = tOl (3.2) 

By Balerkin technique, the vector residue equation for two 
dimensional cutting can be formed as 

JJ Ni |py.W) - 7. g - f I dy. dy = tOl (3.3) 

The term f in equation (3.2) can be included in the pressure 
(isotropic part of stress) by redefining pressure as 

7p ' = S'f f 

The term JJ" Ni (7. dx dy can be written by the weak formulation 
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as 

Expanding the first term in parenthesis by divergence theorem, the 
above equation (3.4) reduces to 

JJ Ni (V. {j;) dx dy = # Ni n.jj; dl - JJ V Ni. dx dy (3.5) 

B 

where the first terms gives the contributions due to the specified 
traction boundary conditions. The original momentum residue 
equation (3.2) can now be rewritten as 


JJ Ni (V. s^) dx 


V (Ni Ni. ^ dx dy (3.4) 


n Ni p (V.V. )y dx dy + JJ 7 Ni. jj; dx dy = # Ni n.j^ dl (3.6) 

--s. 

Expanding g; into its components, the LHS term 7 Ni.jx can be 
written as 


7 Ni.s> 






o 


V 



0.7) 


Momentum equations : 

Putting the above term of equation (3.7) into the 
non jdimensional form of momentum equation (3.6), the x component 
equation can be derived as 

ii2. rf-M- f 

pV t (|Nx fu + V y dx dy 

t ax* ax* y 
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SitTiilarly y -moment um equation can be obtained by taking 

A 

the dot product between unit vector j and vector momentum euqation 
(3.6) . 

In the above equation, the inertial terms have been 
quasi -1 inear ised by using the velocities corresponding to previous 
iteration (or initial guess) values during the iterative solution 
of t|ne problem. 

I 

The field variables have been interpolated within each 
element as shown below 


n 

u = £) Niui 
i = 1 

m 


n 

V = J] Nivi 
i = 1 


n 

T = E NiTi 
I = 1 


/ 


(3.9) 


Similarly, the spatial coordinates of the boundaries of the 

I 

isoparametric element are defined as 


n n 

X = E » y = L Wiyi (3. id 

i = 1 i = 1 

where, n = 1 to 8 
m = 1 to 4. 

It may be noted that the pressure is defined only at corner nodes 
of the element and hence needs only a linear interpolation 
f unct ion. 

All other variables (u,v and T) have been interpolated 
using quadratic shape functions. 

For easy computation, the interpolation function Ni are 
defined in terms of local coordinates as under. 



For corner noes ; 

Ni = ^ (1+ +77.77 -1 ) ( 1+77 -77) 

For mid-side nodes; 

Mi = ^ (1- ?^) (1+77- 77) r = Ox 

Mi = i (i+ ?j?) (1-77^), 77- = 0 i 


C3. 12) 


(3. 13) 


Substituting the ei<pressiQns for the interpolation of 

u,v and p into the momentum equation (3.8) and rearranging, we get 

1 

the form of non -d imensional i sed >;~momentum equation as : 


M 


rii<„ ^ + Rv» "ii ) * ^ 

O' k k ^ k k ^ ^ ^ ^ ^ ^ 


dy. 


dy 


dx dx 


dx^dy^lcu^l 
dy^ J '• 




[rr/t dx^dy^lcp^] = [ 

# Nit^ dl^l 

L-'n ax* ■* J j 1 

J 


(3. 14) 


where a 



a /p\f 
yo 

Traction in 

t fa 
X yo 

Similarly y 


X -d irect ion (non -dimensional ized. ) 


—momentum equation can be obtained 


in the 


form 
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! 


Frr/ .4 ^.4 « 

Cu .1 

J 






' dx’^dy^lcv'fl 



J J 



r# Nity^dl'^j 

(3. 15) 


where and are the traction components per unit area^ given 
by, 

M. 

t = n. O', i 

X ~ ~ 

h 

and t = n. c. j 

y <s» Arf 

The residue minimization principle for the mass balance 
equation takes the form 



dy*-> 


dx*^dy 




O 


C3. 16) 


The weighting functions for the mass balance equation 
have been chosen to be the linear shape functions used to 
interpolate pressure within the element. This is done in order to 
get a pressure field that is compatible with the incompressibility 
condition. Substituting the interpolation expression for the 
velocity variables, the final form of the mass balance equation is 


. [//Mi 


w j ^ ^ w 
ay^ 


dy^j 


Cv^l = 0 
J( 


C3. 17) 


I 


ENERGY EQUATION 
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The weighted residue form of the energy equation is 


J/Ni { 


d , , » \ ^ n ^ . s . w 

- (k ) + (k y } dx dy 

^ V V >C * ' 

Sx dx dy dy 




+ ffNi Pe —— e^dx^dy* 

V3 


JjNiPe (c* u* + c" V* I dx^dy* = 0 (3 

dx dy - J 


18) 


Applying Greene's theorem to 1st term of the above 
equation, the final form of energy residue equation is obtained 
as. 




aNi aNj 


+ k 


ax 


ax 


« aNj^ 

dy* 


aNj 

a * 

dy 


i j+ NiP^ |Cp 


(N uf 
k k 


aNj 


dx 


s 


Kl » 

N, V, 
k k 


aNj 

a * 

dy 


H] 


dx*dy*]cT*J 

= # Ni — — di*+ ^ ~~ e^dx^dy* 13.19) 

an^ ® -Vs 

The second term on the RHS is the heat generation 
contribution and the first term is from heat transfer boundary 
condi t ions. 


3.2.3 BOUNDARY CONDITIONS 

Modelling of boundary conditions is as important as the 
modelling of differential equation themselves. In fact, 
experience reveals that the toughest part of problem formulation 
is the prescription of proper boundary conditions- 


S}^= Convective be. 

Sf = Friction be. 

Sq = Flux b.c. 

1st =Sticking contact 
length Vn = 0 


Sq =0 
, Vc 




U =1 jVsO 
T = Tarr)b 


III t: u» 1 ;v.O 


u =1;v =0 
T = Tamb 

Fig. 3-4 Problem regions showing the boundary conditions 


sticking 


zone 


\ /Normal stress ( 07,*) 


Shear stress (ff*) 


Slipping zone 


Fig. 3-5 Idealized distributions of normal and shear stress 
on the tool rake face 
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The facilitate computation and to implement boundary 
conditions in a systematic manner, the total boundary has been 
divided in to eight separate segments as shown in Fig. (3.4)'. 
Surfaces I, II, and III have prescribed boundary conditions for 
u,v and T. These conditions are implemented by incorporating them 
in nodal equations for corresponding variables, in the final 
matrix equations. Also, the first node of the mesh is prescribed 

with a pressure boundary condition. 

! 

The normal velocity (V = 0) is applied on surfaces IV, 

~n 

V, VI and VIII, since there is no flow in a direction normal to 
these surfaces. Thus, normal velocity at the surface requires. 

V^ = u^ cos 0 i - v^ sin 0 j = 0 ‘ (3.21) 

n 

A A 

where 0 is the angle made by the normal with the x-axis and (i,j) 
are unit vectors in (x,y) direction respect ively. 

Expanding u and v in equation (3.21), 

S * S 

V^ = r N.u^ cos 0 i - E N.v^ cos 9 j =0 (3.22) 

-n j = l J J j=l 

where u^ and V* are the nodal velocities of the element which 
J J 

lies on the concerned boundary. 

On surface VII, which is the far end of the chip, the 
velocity is equal to the chip velocity Vc. The velocity 
components u and v on this boundary are given by 


^ 

tl “ V 

sin a V 


c 

1 

(3.23) 

V = V 

cos o< ^ 


c 
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where, o< = Rake angle 

On the surface V and VI, we apply two boundary 
conditions. 


i) zero normal velocity ie V = 0 

n 

ii) traction boundary condition 


For the normal direction, zero normal velocity condition is 
implemented . 

The traction boundary condition for shear stress can be 
applied by considering the momentum equation in the tangential 
direction as given below : 

(x-momentum equation) (-sina ) + (y momentum equation) ( -cosoi) 

_ fv, . « 

= jNi dj 


where, = shear stress along the surface ( = m/j*e**) 


1 . e. 



^Nj ^ 
ax* 

N. V* 
k k 

^Nj , 
ay 

ax* 

aNj 

dx 

+ fj 

« aNi 

1 — 

ay 

aNj 

ay 

1 dx*dy 

*jsin a. t|U* 

1 

r r r f ^ i 

aNjl 

a;<*J 

dx*dy 

1 

r rFNT'^ 

0< C V . ] 

J 

tTm £«i M.) 
ay." •’ 

dx*dy* 

j sin 

a CP*] 
J 

icc. 

1 

z 

Me 

aNj] 

ay*i 

■ dx*dy*j cos 

« Cu*] 

J 



aNj ^ 
ax* 

N. V* 
k k 

aNj J 
ay* 

, ^ aNi 

+ jj 

a>; 

aNj 

dx 


+2^i* 

aNi 

-a ^ 

dx 

ay*J 

dy cos a 

Cv 1 
J 


;■ 'i 




.V 




j^JJ( M^) dx^dy^j cos o< CPj] - J Ni r^d^ 


(3.24) 
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Tfif? traction boundary condition extended to the 
elements on Vlth surface also, with the assumption that shear 
stress decreases linearly towards the end of the contact patch and 
becomes zero* This based on the observation of many investigators 
(19) that shear stress remains reasonably constant on the sticking 
contact length and reaches zero gradually in the slipping contact 
length* 

For surfaces IV and VIII which have convective heat 
transfer boundary condition, the overall heat transfer coefficient 
is given as input data* 

The weighted residual form of the boundary condition is 




jNik* dl = -J dl = -J Nih"* T^dl + J Nih** dl 

' 


(3.25) 


8 

Expanding T = T] N.T., convective heat lose is given by, 
J = 1 J 


jNik^ — ^ dl = (Nih* N^) dljc T^l - J iMih^ ] 


(3.26) 


On the chip tool interface (surfaces V and VI), the 
condition for frictional heat generation is considered 

This is given by 
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3n 


- = p v:. 
^ e F T 


or 


^n 




?n^ 


3 K, "F 


(3.27) 


weighted residual form of the above equation is 
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|Ni j dl = P ; — V:^dl 
V , - k- 


, The slip velocity is given by. 


(3.28) 


= V cos a + u sin a (3-29) 

I 

ELEMENT ASSEMBLY t 

In the previous section, the elemental contributions to 
the left hand side, coef f icient matriK and right hand side vector 
have been discussed in detail for the eight -noded isoparametric 
quadrilateral elements. These elemental contributions are. 
assembled in to a global matrix equation by adding the entries 
corresponding to each nodal variable appropriately. 

The boundary integral contributions on both left hand 
side and right hand side of the matrix .equation are also 
incorporated as discussed in the previous section. This assembly 
procedure results in a matrix equation of the form 


where 


CAI Cx3 
CA3 
Cxi 
CB3 


CB3 

coefficient matrix 
solution vector 
Right hand side vector 


The problem at this stage is set for solving above 
matrix equation by any standard technique available. 


3.2.4- Matrix Solution Technique : 

For the present problem an iterative technique is used 
for matrix solution based on the Frontal method. Frontal method 
offers many advantages such as small use of memory space and high 
speed of computation. It utilizes the fact that when all the 
contributions for a particular node and for a particular variable 
are over during element assembly, it could be transferred to the 
disk memory. On this basis it retains in its memory only those 
variables which are yet to be assembled the variables are solved 
by Gaussian elimination is after assembly. Since the governing 
equations are highly non-linear an iterative procedure with 
successive under -relaxat ion has been employed, 

3.3 Post processing 

After solving the governing equations, the field 
variables, namely velocities, pressures and temperatures would be 
available at each node of the solution domain. 

These variables can be further processed to get many 
interesting quantities such as forces, stress distribution etc. 
Al4o, it is possible to plot isotherms to study (the temperature 
distribution or iso-strain rate lines which might give an idea 
regarding the size of the deformation rone. 

/ 

Thus, in the present analysis, the following 
post -processing has been done. 

i) plotting of isotherms 

ii) plotting of iso strainrate lines i 

and iii) estimation of forces involved in metal cutting 
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. Isotherms t Temperature has been predicted for all the 
nodes in the domain. Using linear interpolation, it is possible 
to estimate temperature at any point across the element boundary. 
The coordinates of points of with same temperature are calculated 
using this procedure and joined using smooth curves. At a simple 
glance, isotherms give a clear idea about the temperature 
distribution in the domain. 

. Isostrainrat© Curves : Once the velocity fields all over 
the domain are available, the strain rates can be evaluated veri^ 
easily using equation (2.14). The strainrate plots give a good 
idea about the deformation pattern in the domain. The strain 
rates at different points are first normalized with respect to the 

I ' 

maximum strain rate occurring in the solution domain. Points 
haying the same value of normalized strain rates are connected by 
smooth curves and the curves 2X, 3X, 4X, 5X, 6X, 8X etc. are thus 
plotted. These constant strain rates provide a clear picture of 
primary and secondary deformation zones. 

. Estimation of Forces * Forces acting on .the tool are of 
paramount importance. Now, once the velocity fields, pressure 
fields etc. are available in the domain, finding stresses and 
forces is not difficult. For force calculation, only the elements 
lying on the rake surface would come into the picture. (i.e. 
surface V and VI). For there elements the stresses can be easily 
calculated using the following relations : 


56 


a 


Y 

XX 


-P* + 

dx' 




a 


yy 


-p + 2 m 

dy/ 




. ^ W 

and <y = o' 

xy yx 


£h_ + 


Further, these stresses can be used to find the shear 
stress and normal stress on the rake surface, using. 


,a . O' N /O 

a ( X X + y y ) ( x : : 

nr) = ^ Lr,. — +. 


a 


Ll — L cos2o( + sin2o( 


xy 


and 


» « 

,o O' . 

» -( XX - yy ) . „ . » _ 

T_ = n — sin2o( + O' cos2o( 

F 2 xy 


w; W » “St 

Tp also can be calculated using the relation, Tp = m ffs c 
where, m = friction factor 

e = stranrate teffective] 

■= Normal stress acting on the rake surface 
nn 

- shear stress acting on the rake surface. 

F 

From these, we can easily calculate normal force and 
shear force acting on the rake surface. 


F ■ = 


n 


where , 


n 


X contact length x width 
=\r^ X contact length x width. 
= normal force 


57 


F = shear force 
s 


and F^ can be easily resolved to gel the desired 


forces viz F^ and F^ 


F^ = 

X 

- f’^ 

S 

si nr 

- F^ 
n 

cosr 

F^ = 
y 

_ 

n 

sin?'' 

- F^ 
s 

cosr 


where, F^ = cutting force or_ thrust force 


- feed force, 

y 


y = rake angle 

These forces can be expressed in dimensional form as. 


F = f"** >; O' X t^ 

X X yo 


V 2 

F = F ■ X O' X t 
y y yo 


wherer <:>* = uniaxial yield stress of the workpiece material at 

yo 

ambient temperature 
t - Depth of cut. 

SincCy we are using iterative solution technique, for the 
calculation of forces, the fully converged values of the variables 
are used. 


3. 4 Programming 

Computer programming is one of t\\B challenging tasks of 
the present work. 

Tite pr agr ammi rtg work is divided in two parts. 
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(i) program for grid generation. 

(ii) program for finite element analysis 

' I 

3.4.1 Program for grid gonoration t 

' A simple flowchart shown in Fig. (3.6) gives a general 

picture of the procedure involved. 

In the main programme itself, the data related to the 
cutting process such as depth of cut, rake angle, contact length 
and chip thickness are read. If contact length and chip thickness 

are not available, provision has been kept to calculate them 

within the program. All the distances have been 

non -dimensional i zed using the depth of cut. Other relevant data 
for the domain such as number elements in >: direction, and 
/-direction tfor each zone), number of common elements, number 
nodes per element etc. is also read in the MAIN itself. 

Subroutine ZONE is called individually for each zone to 
generate nodal connectivity and them a globalised connectivity is 
generated using subroutine SUPER after considering common 
elements. 

Subroutine AUTO is called to generate co-ordinate data 

for each node, separately for each zone and later this data is 

global ised . 

Coordinates of the corner points of the zones and of 
some specific points on the boundary Ceg. the point from which a 
finer mesh is desired etc.) is supplied in subroutine INPUT. 
Using this data, first, coordinates for all the boundary nodes is 

' I 

generated. Generation of these coordinates is simple so long as 


the boundaries are straight. 


For the eighth, sur face, where the concept of hyperbolic 
stream line is used (explained elsewhere). Subroutine CURVE is 
called, to generate the boundary nodal coordinates along this 
curved surface. 

All this coordinate data for the nodes on the 
boundaries, is used to generate the coordinates for all the 
elements in the domain in the subroutine INTPOL, where the concept 
of trans-finile interpolation (explained elsewhere) is used. This 
is done for each zone separately. Finally, all these datas of 
individual zones are globalized. 

3.4.2 F’rogram for Finite Element Analysis 

Fig. 3.7 gives the flow-chart of the program. 

' I 

Main program calls the subroutine DIMENS, DINPUT, HERAT 
and FORCE, which inturn call the other subroutines. 

/ 

Subroutine DIMENS returns the values of the array 
dimensions which remain unchanged throughout the program. 

Subroutine DINPUT reads and returns all the input data 
such as cutting conditions, initial conditions, boundary 
conditions, nodal coordinate data, element connectivity data and 
prescribed parameters like ambient temperature and material 
properties at ambient temperature. The whole input data is in the 
non-dimensional ized form. Some of the inputs are checked for 
correction by two routines named DIAGMl and DIAGN2. 

Subroutine HERAT calls the FRONTS routine for solving 
assembled matrix equation. This routirse begins the problem with 



Fig. 3-7 Flowchart of Finite element analysis program 


z o > 























guessed initial values for all the Tiodal variables and relaxes 
1 

them for next iteration in the subroutine TOLREli, if the solution 
is not found converged within the prescribed tolerance. This 

routine also calls the WRITER subroutine (i) to write the nodal 

/ 

values of the variables if the solution has converged and Cii) to 
write the largest relative change in the variables that has 
occurred during a particular iteration if the solution is not 
converged. ; , 

Assembling is MATRIX subroutine which calculates the 

I 

element fluid matrix flumx. This fluid matrix consists of the 
elemental contributions of the governing equations to form 
coefficient matrix. Similarly elemental right hand side vector is 
also calculated here. Subroutine DRIVES is called within MATRIX, 
which in turn calls the subroutines SHAPES, DJACOB S: SHAPE4 to 
calculate the shape function and their derivatives. 

The boundary conditions are implemented in BOUCQN which 
is called in MATRIX before it returns to FRONTS., This is called 
for only boundary elements. Here, relevant boundary conditions 
from among the zero normal velocity, convective heat loss, the 
frictional heat flux of frictional shear stresses are identified 
appropriately and elemental contributions are evaluated. The 
contributions to right hand-side vector are calculated at local 
nodes and are assembled in the global vector both in MATRIX 
BOUCQN subroutines. 

Subroutine FORCE does the post processing part of the 
nodal values predicted after convergence. Frictional force and 
the normal force on the rake surface are evaluated considering the 
elements lying on this surface and finally, the cutting force and 
the feed force are evaluated from them. 


CHAPTER 4 


EXPERIMENTAL ANALYSIS 


When we speak of theoretically estimating the quantities 
of practical importance such as forces and temperatures, it is 
obviously needed to check these predictions with those obtained 
experimental ly. 


The following experiments have been conducted in the 
present work. 


Measurement of average temperature along the chip-tool 
inter f ace. 

Measurement of Cutting Forces. 


The tests have been designed for various cutting 
conditions to find the effect of cutting speed, rake, 
feed etc. on the temperature distribution and Forces. 


4.1 TEMPERATURE MEASUREMENT 

I , 

The determination of temperature distribution on the 
rake face and in the cutting tool is very important for estimation 
of tool wear and tool life. Several techniques have been used by 
researchers for accurate determination of temperature at different 
points on the tool. Some of these are 


(i) Use of photosensitive paints, 

tii) Use of tool workpiece thermocouple. 

(iii) Using tiny thermocouples which are embedded in the tool 

and 

(iv) Using infrared photographic technique. 
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Most of these techniques except the tool -work 
ther mocoup 1 e technique are not popular because of the low spread 
of response or complexity. Although tool -wor kpiece thermocouple 
procedure gives only the average temperature, it is a very simple 
and reliable method. This method has been employed here. 

4. 1 • 1 Tool— Wor kpioc© Thormocouplo Method 

\ 

I 

This method makes use of the emf produced between the 
hot contact patch of tool and workpiece and their cold ends. If 
the cold ends of the tool and workpiece are joined, a smail 
current will flow which can be measured. Also, the emf generation 
can be measured in terms of volts. 

This technique is based on the following laws of thermo- 
electric circuits (MCS). 

I 

1. The emf of thermoelectric circuit depends only on the 

difference in temperature between the hot and cold 
junctions, and is independent of the gradients in the 
parts making up the system. 

2w The emf generated is independent of the size and 

resistance of the conductors. 

3. If the junction of the two metals is at uniform 

temperature, the emf generated is not affected, if a 
third metal, which is also at the same temperature, is 
used to make the junction between first two. 

A great number of thermocouple hot' junct ions are active 
at the chip-tool interface asperities which are forming a parallel 
electric network. Thus, the emf recorded would be an algebraic 


65 


average of all the thermal emfs. This average emf can be 
calibrated against temper ature and the average temperature can be 
found by this technique. 

4.1.2 Tost Conditions 

All the experiments have been carried out on a HMT 
lathe. The workpiece material selected was Mild steel. HSS tools 
of standard dimensions with following rake angles have been used. 

lO"", 15®, 20® 

The tool was held firmly in the position of tool post 
within the tool holding arrangement provided on the dynamometer. 
A constant overhang of the tool of 18 mm was maintained for all 
the tests. The experiments were carried out under dry cutting 
condit ions. 

Due to some practical limitations of machine and the 
measuring devices, the experiments were conducted only for lower 
values of feed (e.g. capacity of loadcell —* 100 Kg.) <,0.05 mm/rev 
to 0.075 mm/rev). The cutting speed was conveniently selected in 
the range 0.4 m/s to 0.8 m/s. 

Rake angles of 10®, 15® and 20® were used. The 
fO|Howing cutting conditions were maintained constant for all the 
tests : width of cut = 2.4 mm and clearance angle = 10®. 

However, only a few of these readings taking care of the 

I 

variations in all process parameters viz. feed, speed and rake 
etc. have been used for comparison. 
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4.1.3 Experimental Procedure 

The block diagram of experimental eet-up is shown in 
Fig. (4.1). The workpiece is insulated from the chuck using a 
insulating sheet. Even the tool holding arrangement within the 
dynamometer is insulated perfectly. Care has been taken that 
contact exist only between the tool and workpiece at the cutting 
edge. The workpiece and the tool tip constitute the hot junction 
H. The cold junctions A and B remain at room temperature. 
Mercury is used to make electrical contact with the rotating 
member in order to complete thermo-electric circuit. The wires C 
and D may be ordinary copper wires since both the ends are at the 
same temperature and these have no influence on the emf generated. 
The emf generated is recorded by a double channel recorder/ 
plotter • At the end of each experiment the tool and the workpiece 
were allowed to cool to the room temperature. This was done to 
ensure that the tool tip is always at room temperature at the 
beginning of the experiment. In each of the test cases, the emf 
was recorded after a stabilized trend was observed. 

The measured values of emf were converted into 
temperature with the help of a calibration curve. The calibration 
wa's carried out by contacting the tool tip to aiMS plate put in a 
furnace with sufficient insulation. The temperature of the hot 
junction was measured by a standard Chromel -Alumel thermocouple, 
which was later used for comparing with corresponding emf readings 
of tool- work thermocouple. Thus the calibration curve was 
prepared by varying the temperature of the MS plate in the 
furnace. 



/ 



Fig. 41 Exptl. setup for the measurement of average tool-chip 
interface temperature and forces. * 


Fig. 4-2 Schematic view of dynamometer for measuring cutting 
forces in two dimensional turning operation. 
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4. a FORCE MEASUREMENT 

4. a. 1 Rolationships between various force coiT 5 >onents s 

The notations used for the forces are as follows. 

i> Main cutting force coinciding with the principal cutting 

velocity direction is taken parallel to X-axis of the 
Cartesian coordinate system and is denoted as F (Fig. 

' 4.2). ' 

ii) Feed force acting in perpendicular direction to main 

cutting force in Y-direction, is denoted as F^ (Fig^ 
4.2). 

The action of forces in a single point orthogonal two- 
dimensional cutt ing operation producing cont inuouschips can bp 
readily understood from Merchant's circle diagram shown in Fig. 
(4'. 3). Here F^ and F^ are Cartesian force components acting upon 
the tool nose, normally called main or tangential component and 
radial or thrust component. The shear plane of the material 
sustains the normal and the shear component (F^l_^)- For 

force balance, from free body diagram of the chip we get the 
conditions that the resultant force (R) on the rake surface should, 
be CO-1 inear and equal in magnitude to that acting over the shear 
plane. Thus the circle with R as the diameter will contain all 
the force considered, and this diagram can be used conveniently 
for the mutual conversion of the component forces into equivalent 
shearplane and rake forces. Thus we can write. 
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^h 

r 

cos<^ 

- F 

y 

sin<;6 

^h 

= F 

y 

cos^ 

+ F 


F 

s 

= F 

X 

sin^^ 

+ F 

y 

COS^' 

F 

n 

= F 

X 

cosy- 

- F 

y 

sin^ 


If the friction coefficient on the rake face is fJr then 

we have 

= tan ^ fj 



n 


where I’i is the friction angle. 


In practice the external cartesian force components 
and Fy are usually measured which is then made use of in further 
force evaluations using the equations already derived. 

It is well known that force is an entity generally 
measured through its action on a system offering a finite 
resistance, which is usually pre-calibrated. Based upon the 
principles of modern transducer techniques, the measurement of 
cutting force hence involves, a transducer, signal conditioner 
which converts the transducer output into a useful measurable 
parameter and lastly an indicator or a recorder to register the 
f i n a 1 e s u 1 1 s . 

Tfie accurate measurement of F and F (components of 

1 ^ , y 

resultant tool force) has been the subject of considerable effort 
in the past, and several types of cutting force dynamometers have 
been develcsped. In most of the metal cutting force dynamometej^s 
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the tool force is determined by measuring the deflections or 
strains in the elements supporting the cutting tool. It is 
essential that the instrument should have high rigidity and high 
natural frequencies so that the dimensional accuracy of cutting 
operation is maintained and tendency for chatter or vibrations to 
occur during cutting is minimized. They should have reliability 
over a range of temperatures. 

4r. 2. lost- Conditions for Force Measurements 

Test conditions were maintained same as those used in 
temperature rrieasur ement . 

I 

The whole set of readings namely, 

i) Force F 

K 

ii) Force F / 

y 

iii) Temperature 

iv) Chip thickness 

were taken for various combinations of cutting conditions. 

The dynamometer used in the experiment was cantilever 
l!ype straingauge mounted. Concept of Wheatstone bridge was used 
in the circuit. The constructional details of the dynamometer is 
not being explained herer since it is easily available in any book 
on metal cutting. The capacity of the dynamometer was 100 Kg. 

Omniscribe stripchart recorders were used for recording 
forces as well as temperature readings. 

4.2.3 EXPERIMENTAL PROCEDURE 

The experimental set up for orthogonal cutting is shown 
in Fig, (4.1). The workpiece selected was a mild steel seamless 


f 


tube of suitable dimensions. To begin with, the ID and OD of the 
tube were smoothened by taking some cuts, till thickness was same 
throughout . 

After setting particular combination of cutting 
conditions, forces and temperature readings were taken 
simultaneously and some chips were collected. Chip thickness was 
found by using the weight method. For some of the readings, 
attempts were made to measure the contact length by applying lamp 
black on the rake surface and then, measuring length of scratch 
under a travelling microscope. 

Many researchers have done a lot of experimentation of 
this type. We have used some of this work for comparing our 
results. Temperature distribution in workpiece and chip during 
orthogonal cutting was observed by Boothroyd 1191 using infrared 
photograph technique. This has been used to compare the isotherm 
pattern predicted from our analysis. Muraka et al. C173 have also 
done a compr ehensive experimental work consisting of measurement 
of forces in metal cutting. This has been used to compare the 
predicted forces from our analysis. 

A summary of experimental results used for comparison is 
given in Table 1- 




f 
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'..I ’ t, r •* q r f-i n H -5, i, n v] " • n fi M n r\ s \.{ r r-: d D a t a 


No. 

W'')r 1: niece 
Revs/rriiri 

N rom 

[ p ifJr c 

diartipler 

D iwn 

Width 
of cut 
W fTin 

'Speed 

V 

rri/rTiin 

Feed 

t 

tiiro/ 

rev 

Ra^-.e 

Angle 

iQ) 

Chip th; 
nesB rai 

r 

icF- Shear 
:io angle 

0 

(0) 

Contact 

length 

>t 

m\ 

Cutting 

force 

Fx 

N 

Feed 

Force 

Fy 

N 

Average Ct| 
tool inter'! 
face terip.;r 
°K 1 

1 

57 

162.10 

2.35 

29,02 

0.05 

10“' 

0. "^9 

27. B 

0.14 

563 

3S9 

479 j 

n 

£ 

57 

167.10 

2.35 

29.02 

0.063 

10° 

0.434 

27.4 

0.17 

641 

416 

484 

3 

57 

162.10 

2.35 

29.02 

0.075 

10° 

0.5 

29.3 

0.2 

702 

452 

502 j 

4 

71 

162.10 

2.35 

36.15 

0.075 

10° 

0.49 

27.3 

0.2 

667 

453 

II 

500 ii 

''1 

5 

?0 

162.10 

2.35 

45.83 

0.075 

10° 

0.46 

26.2 

0.2 

633 

449 

502 

+ 6 

80 

98.20 

6.35 

24,68 

0.3556 

20° 

0.536 

39° 

1.11 

3525 

1560 

- 

+ 7 

80 

98.55 

6.35 

24.77 

0.3556 

30° 

0.61? 

qcO 

0.7366 

2650 

545 



Note : 

* All the tests carried out in dry cutting conditions, 

, For all the tests HSS single point cutting tool is used, 

BticHng contact length, 

E;:per ifTjsntal results of fluraka et al. C17]. 


^ I 
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CHAPTER 5 


RESULTS AMD DISCUSSIONS 

the present work, a finite element formulation for 
analysis of orthogonal metal cutting has been 

Based on this formulation, a complete software 
package, including useful post and pre-processing modules has been 
provided. 

The prirriary objectives of the present analysis have been 
as follows : 

i) To predict velocity, pressure and temperature fields in 

vicinity of tool tip, by the finite element solution of 

^ coupled stress balance and heat balance equations. 

■ I 

ii) To conduct a detailed parametric study, to highlight the 

effects of process parameters such as, cutting velocity, chip 

thickness ratio, depth of cut etc. > 

iii) To predict quantities of practical importance such as forces 
involved in metal cutting, isothermal patterns and the shapes 
of deformation zones. 

I • J 

iv) To validate the theoretical model with experimental results. 

! 

v) To verify the application of energy minimization principle 
for determining the value of chip-thickness ratio. 

In the course of analysis, some significant quantities 
such as strain rates, which are very important for the study of 
the deformation pattern in the plastic region, emerge out. The 
variation of strain rate distribution across the shear plane and 
along the shear plane have been studied to throw light, upon the 


In 

the thermal 
developed . 
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ds f ov iTiri t i on ptocsssss undergons by thn iriaterial. Also, ths 
iGmpHrciturp distributions siong ths rsks 'f‘ 3 cs in ths con 1 3 C t 
rsgion and in ths vicinity of ths shsar plane have been plotted to 
analyse the thermal phenomena in detail. 

5.1 Isotherms and Deformation Patterns 

The temperature values at all nodes have been predicted 
through the solution of the energy balance equation. Isotherms 
are then plotted by interpolating the available nodal 
temper at ur es r linearly on element boundaries, to obtain the 
locations where the given temperature values (isotherms) occur. 

The strain rates (e) at different points are first 
normalized, with respect to the maximum strain rate occurring in 
the solution domain. Points having the same value of normalized 
strain rates (after suitable interpolation) were connected by 
smooth curves and the curves of 37., 47., 57; and BV. normalized 

I 

I 

strain rates (as compared to the maximum of 1007.) were thus 
plotted. These constant strain rate curves give an idea about the 
size and shape of the deformation zone. 

I 

Isotherms and deformation patterns for various cutting 
conditions are shown in Figures 5.1 to 5.5, 

It is observed from the isothermal patterns that, as 
expected, the maximum temperatures occur near the chip-tool 
contact region. This has been observed by earlier researches as 
well. The reason for this is that, intense heating of a material 
by frictional work occurs at the contact region, in addition to 
earlier heating due to plastic work in PSDZ. It is also seen that 
temperatures in PSDZ are quite less than those occurring in SSDZf 




Temp in®K 



Fig. 5-2 (a) Isothermal pattern (b) Deformation pattern 





Fig. 5-3 (a) Isothermal pattern (b) Deformation pattern 
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from this trend, it is clear that in metal ' cutting problems, 
frictional heat generation has an important effect on the 
temperatures. Bulk of the material, especially the domain nearer 

I 

to uncut material is at room temperature. This can be attributed 
to the fact that a major portion of the heat produced is carried 
away by the chip. High temperatures are also observed in the 
freshly machined p-^rt of the workpiece, just below the tool tip. 
Some of the heat generated in PSDZ and SSDZ is conducted to the 
uheut material causing high temperature values just below the 
cutting edge; moreover, the relaxation processes which occur in 
that region also result in some heat generation. The occurrence 
of very high temperatures on the rake face contributes a lot to 
the tool wear. 

Deformation patterns obtained in the present analysis 
confirm the well-established observations about the existence of a 
narrow PSDZ bounded by two parallel planes and a wedge-shaped SSDZ 
adjacent to the rake face. Among the various iso-strain rate 
curves, it appears that 57 . cut-off gives a reasonable 
correspondence with the PSDZ thickness obtained through an 
established relation proposed by Jain 8c Gupta £27). 

It is observed that the constant strain rate curves on 
the uncut material side are closely spaced compared to the cut 
material (i.e. chip). The reason for this may be that the 
material approaching the tool tip undergoes deformation near the 
shear plane, due to the shearing action caused by high compressive 
forces of the tool. For the chip the deformation occurs over a 
wider region since one of the surfaces is unrestrained. In a very 
small region of the uncut material i.e. just below the tool tip, a 


close net of e>;tende'j iso-strain rate lines is found. Since 
separation frorri the chip material occurs in this region, high 
stress values may be prevail which explain the above trend. The 
secondary deformation zone appears to be of triangular shape again 
confirming, a well established observation. Close -spacing of 
strain-rate lines in SSDZ is again attributed to high stresses in 
that region. 

S. 2 Effect of Process Variables 
5. 2. 1 On tool Face ten^eratures t 

It is observed from Figs. 5.6 to 5,8 that a non-uniform 
temperature distribution exists on the rake face and the maximum 
temperature always occurs at some distance away from the cutting 
edge. It can also be seen that, the temperature on the rake face 
rises sharply up to a point of maximum temperature with distance 
measured from the cutting edge and then drops moderately. Near 
the end of contact patch, the temperature more or less remains 
constant. This pattern has not been observed by earlier 
researchers. Muraka (17) observed that the temperature drops from 
a rrnaximtim near tool tip to a minimum at the point of chip 
separation. He had also observed a sharp change in the 
temperature gradient. Sarma (21) observed the occurrence of two 
local maxima. These differences may be attributed to the ignoring' 
of effects such as thermal softening of material variations in 
thermo physical properties and friction factor variation along the 
contact patch, in the previous works. 

An investigation into the effect of rake angle on rake 
face t etTiper a t ur e reveals interesting trends. With increase in 
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Fig. 5-6 Effect of depth of cut on rake face temperatures 
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Fig. 5 -7 Eff 





Temperature PK) 



Fig. 5-8 Effect of rake angle on rake face temperatures 



down. 


This is 


rake angle, the temperature distribution falls 
because, the included angle (i.e. angle between the tool flank and 
face) reduces with increase in rake angle. In turn, the total 
energy input to the ^ystem is reduced and hence the amount of heat 
generated is also decreased, thus causing a reduction in 
temperature. 

Increasing the depth of cut or velocity causes an 
increase in the temperature distribution, proportionally. The 
reason for this is that increase in depth of cut or velocity leads 
to an increase in the energy input to the process. 

5. 2. S On Temqporatur© distribution in Shoar Zon© 

Figures 5.9 to 5.12 show the effect of various process 
parameters on temperature distribution along and through the shear 
plane. The trends with respect to cutting parameters seems to be 
similar to those seen for the rake face temperature. The 
temperature values in the region near the tool tip are much higher 
than those near the free end. This can be attributed to the 
decrease in thermal conductivity of the material with increase in 
temperature, which restricts the outward flow of the heat 
generated in BSDZ for PSDZ. Secondly, the free 'surface being in 
contact with the atmosphere convective heat loss is expected. In 
the present work, shearplane is established as parallel to 107. 

I 

cut-off strainrate curves, and in the mid of PSDZ, approximately. 

S. 2. 3 Variation of Strain-rates in Shear Zone 

Strain ratp variation along and across the shear plane, 
have been drawn in Figures (5.13 and 5.143. As can be expected, 
acr'oss the stiear plane strain —rate rises rapidly, reaches a 
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Fig. 510 Effect of cutting velocity on temperatures along 
shear plane 
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Fig. 5-14 Strain rate variation normal at shear plane 



maxinfiuiTi somewhere near the shear plane and falls rapidly again- 

! 

This confirms the ejtistcnce of a narrow 5SDZ near the shear plane. 
Also, it is observed that the strain rate is very high near the 
tool tip and less near the free surface o'" the middle part of the 
chip. Such a trend is to be anticipated, since SSDZ is also 
present near the tool tip. 

Thus the results of present analysis have confirmed all 
the well-l^nown trends in metal cutting to a reasonable extent. 
The width of PSDZ showed a slight decrease with the increase in 
cutting velocity, as observed in experimental studies. 

5.3 Con^arison with experimental Results 

The primary objectives of the present analysis have been 
the prediction of cutting forces and the validation of the 
theoretical model using experimental results. 

A comparison of the forces and temperatures obtained 
exper imental 1 y and .theoretically is shown in Table. 2. 

The forces (both cutting and feed forces) and 

temperatures have been predicted reasonably well. The analytical 
results deviate from experimental values just by 18-2551. 
Considering the complexities involved in the problem, the 
prescription of boundary conditions and the inherent, limitations 
of |:he visco-plastic model, the agreement seems to be quite 

satisfactory. The predicted temperature values are some what on 

the lower side when compared with the experimental values. This 

may be due to not having precise values of the friction factor on' 
the chip-tool interface and heat transfer coefficient on the chip 
surface. It should also be kept in mind, however, that the tool 


92 


Table 2, t Comparison of analytical and experimental results 


Test 

No. 

Cut ting 

(F > 
y, 

N 

Force 

Feed Force 

CF ) 
y 

N 

Average chip-tool 

interface temp. 


Exper i - 
mental 

Anal y - 
t ical 

Experi- Analy- 
mental tical 

Exper i - 
mental 

Analy - 
t ical 

1 

563 

461 

389 

313 

469 

398 

2 

641 

582 

416 

346 

481 

412 

3 

702 

562 

452 

351 

502 

445 

4 

667 

535 

453 

349 

500 

456 

5 

633 

514 

449 

334 

502 

458 

6 

3525 

3618 

1560 

1820 

- 

463 

7 

2650 

2258 

545 

869 

- 

561 


1 


I 


/ 



1 



Total work done (x 10 W) 


Test no. 6 t = 0- 3556 i Rake = 20° jV = 24-68m/min 
Chip thickness 
Exptl. measured =0-67mm 

Merchant’sl \ 

f 0-637mm \ 

J \ 

Present =0-62mm 

analysis — 


C hip thickness (mm ) 


Fig. 5-15 Prediction of chip thickness from energy minimization principle 
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work thermocouple method used for experimentation does not give 
accurate results. Fig. 5.0 shows that ttie isotherms predicted 
reasonably match with the those obtained by Boothroyd tl9) using 
infrared photographic technique. Slight deviation may be due to 
the above mentioned reasons. 

5. 4 Energy minimization Principle 

The Classic Merchant's theory uses the energy 
minimization principle in establishing shear angle or the chip 
thickness ratio. 

I 

I 

In the present analysis, the total work comprising of 
total plastic work and frictional work is obtained for various 
chip thicknesses. This was done generating separate grids for, 


various chip 

thickness 

values. 

f ol lowed 

by the 

finite element 

analysis. It 

has been 

observed 

that 

the 

total 

work falls 

down 

with increase 

in chip 

thickness 

upto 

a 

certain 

value and 

then 

again rises. 

! 

Thus, it 

can be stated 

that there 

exists a 

chip 


thickness value for which the total work or energy involved is a 

I 

minimum. This chip thickness value matches reasonably with the 
predictions of the Merchant. The variation of total work with 
chip thickness is plotted in Fig. 5.15. 

Thus, the present work provides an experimentally 
validated theoretical formulation and software package for the 
visco-plastic and thermal analysis of orthogonal metal cutting. 


I 

CHAPTFR 6 

CONCLUSIONS AND SUGGESTIONS 

From the detailed finite element analysis of the orthogonal 
metal cutting^ the following conclusions can be derived. 

I 

i) Viscoplastic model of the metal cutting can be successfully 
analysed using FEM to study the deformations and thermal 
phenomena thoroughly. 

ii) Borne significant quantities of practical importance such as 
forces, temperature distributions etc. can be predicted, 
accurately. In fact, hardly any work has been done 

especially for estimating forces in metal cutting. 

iii) The constant strainrate curves plotted with cut off values 
of 57-, 47-, 37. and 2.7. give a reasonably good idea about the 
sizes and shapes of the deformation zones. 57. cut off gives 
a better of PSDZ thickness. 

iv) The effects of thermal softening and variations in material 
properties with the temperature are very important. After 
considering these variations only, the isotherm plots 
started showing good agreement with experimental 
measurements. 

v) The parametric study conducted shows that the effects of 
depth of cut, cutting velocity and rake angles on 
temperature of strainrate distribution are significant. 

vi) The thermal model developed confirms the energy minimization 

' principle of metal cutting as proposed by Merchant for the 
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determination of the chip thickness ratio and shea^ angle. 

Suggostions for Futuro Work 

i) Similar analysis may be extended for a detailed stress 
analysis of the workpiece and tool, 

ii) Tool wear phenomena can be analysed thoroughly using the 
temperature and force predictions. 

iii) Effect of contact length on the deformation phenomena can be 
studied by using some of the new concepts such as controlled 
contact . 


I 
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D I MEM3 1 ON COORD ( 386 , 2 ) 

COMMON/Al /LNODSC 111,8) ,LN0DB1 (27,8), LN0DS2(84, 8) 
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FI =3. 142/4+0. 5»(RAKE-BEFA) 

CHR=-TAN(FI ) / (CnBA+SINAwTAN(FI ) ) 

CHPTHK=DCUT/CHR 

CDMLEN=T)CUTr*SI I N ( F I +BET A -R AK E ) / ( S I N ( F I ) sCOS ( BETA ) ) 

WR I I E ( St , 9 ) CHR , CHPTHK , CONLEN 

9 FORMATT.//' CALCULATED VALUES ARE ' 

! //' CHIP THK RATIO ',F8.4, 

I /' CHIP THK = ',F8.4, 

! /' CONTACT LENGTH = ',FS.4) 

PAUSE 'PRESS RETURN TO EXECUTE' 

GOTO 555 
444 WRIT F'^:, 3) 

3 FORMAT (/' ENTER CHIP THICKNESS') 

READF.', , :>) CHPTHK 
WRTTE(;'t,4) 

4 FOI/MATT,/' ENTER CONTACT l.,ENGTH - ) 
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r^,FAD<.w, ?^,K:;(jrii..EN 
555 Cl IPT 1-lK-CHF^THK/DC'JT 
COMLFFN^CONLEN/'DCUT 
MX 1^9 
MY1"=3 
MX2^14 
MY2-^'6 
MN0Dr’-:8 
MEI..Cf:,]f'1--6 

N[.)l''C:uri^--"2-:"4IEl..C0M4 1 

Nl-:LEMC=NXlf^NYl 

ME!..,E,M2^^NX2'«MY2 

MEL Ert-NELEM 1 . +NEL.EM2 

MP C) I I'l 1 =• (, 2«N X J, + 1 ) w ( 2^=;NY 1 + 1 ) -N X 1 ?-MY 1 

MP0IM2™(2wNX2+l) w(2wNY2'-l ) -MX2wNY2 

NPOIM-^NPOIMI H4POIN2-(NELCQM«2+l) 

DO 300 I=n,NPDIN 
COORD (.1 , 1 )=-0,0 
C00RDa,2.U--0.0 
300 COMTIMUE 
ICALL~1 

NT0TC1“2?;NPDIN1 
NBC0M1"(NX1 fNYl)w4 

CALL AUTO (NX 1 , IMYl , NBCONl , XNCIDl , YMODl , NTOTCl , I CALL) 

ICALL^ICALLM 

NT0TC2"'2«NP0IN2 

NBCQN2== (NX2+NY2 .) w4 

CALL AUT0(NX2, NY2, NBC0N2, XN0D2, YN0D2, NT0TC2, I CALL) 
CALL SUPER(LM0DS,NSIF1) 

DO 1.0 I = 1,NP0IN1 
COORDdr 1)=XN0D1(I) 

C00RD(I,2)=YN0D1(I) 

10 CONTINUE 
JJ^-0 

DO 20 I=N0DC0M+1,NP0IN2 
JJ^==JJ-fl 

COOFID ( MPO I N 1 -d J , i. ) =XN0D2 ( I ) 

COORD ( NP 0 1 M 1 + JvT , 2 ) =YNaD2 ( I ) 

NrTMF’DIML+ JJ 
20 CONTINUE 

WRITE (29, 56) 

56 FOF'MAT (// , 6X , "'MnDE-' , 3X, ■'X-COORD-' , 4X, ■-Y-COORD-' , / > 

DO 1239 L-=1,NP0IN 

COORD ( T , 1 ) :: COORD (1,1) ^^DCUT 
COORD ( 1 , 2 ) '-TOJORD ( 1,2) f^DCU I 
WR I T r ( 29 , 1 08 ) I , COORD (1,1), COORD ( 1 , 2 ) 

1 03 f ■ ORM A r ( 6 X , 1 5 , 2F 1 0 . 4 ) 

1289 COM'IIMUE 

TYPE.-:i, ■PL.EIEF.-: COORD. OUT FOR COORDINATES;-' 

STUf'-' 

END 

C 

SUBROU'r IME SUPER(LNODS, MSIP 1 ) 


c 


c 


c 


D I MF: MB I ON LMODB (111,8), LMODS 1(27,8), LN0DB2 ( 84 , 8 ) , LNODSX (111,8) 
COMFKJM/ A 1 1 /NX 1 , NY 1 , NELEM 1 , NPQI N 1 
CGNMON / A 1 2 / N X 2 , NY2 , NELEM2 , NP □ I N2 


DO 600 I =1, NELEM 1 
DO 600 J=l,8 
600 LNDDSl (I, J)=^0 

DO 610 1=1, NELEM2 
DO 610 J=l,8 
610 LN0DS2(I, J)=0 

DO 630 1=1, NELEM 
DO 630 J=l,8 
630 LM0D3(I,J)=0 

WRIIE(24,26) MELCOM 

26 FORMAT ( / 1 0 X , " NO . OF COMMON ELEMENT S = ■- , 1 4 , / ) 
NP0INX=NP0IN1 
I NELEMX--~-NFLEMl 

CAL. L Z ONE (. NX 1 , NY 1 , NPO I NX , NELEMX , LNODSX ) ' 

DO 99 KK=1, NELEMX 
DO 88 JJ=l,a 

LNODS 1 (KK , JvT) =LNODSX (KK , JJ) 

88 CONTINUE 
99 CONTINUE 

WR I TE ( 24 , 22 ) NX 1 , NY 1 , NELEM 1 , NPO I N 1 
22 FORMAT (//5X, -'NXl^', 13, 2X, 'NYl^', 13, 2X, ■'NEI.EM1=' , 15, 
!2X, ■••MP0IN1 = 'M5, /) 

WRITE (24, 33) 

33 FORM AT (//I OX, "ELEMENT CONNECTIVITY INFO. FOR ZONEl", 


! //,6X, ■■ lELEM", 16X, -'LOCAL NODE NOS",/, 

! 16X, ■"1-',5X, "C'^SX, •'3-',5X, •'4',5X, •'5-',4X, "6",6X 
DO 500 IELEM=1,NELEM1 

WRITE (24, 11) lELEM, (LNODBl ( lELEM, JJ ) , JJ=1 , 8) 
11 F0RNAT(5X, I5,2X,8(X, 15)) 

500 CONTINUE 


•' , 5 X 


NP0INX=NP0IN2 


NELEM X=NELEM2 


/ 


"S', / .) 


CALL ZONE (.NX2, NY2, NFOINX , NELEMX , LNODSX ) 


DO 999 KKK=1, NELEMX 
DO sea JJJ=1,3 

LM0DS2 (KKK , JJJ)=LNODSX (KKK , JJJ) 
,838 CONTINUE 


999 COMT ]' NUE 

WP.T 1 F (24 , 23) MX2, NY2, MELEM2, Nr’01l-I2 
23 F0RMAT(//5X, -"NXS^", 13, 2X, -"NY-Z^", 13, 2X, 'NELENZ-:" , 15, 

!2X, "NP01N2Y- ' , 15/) 

WRI rE(24, 34) 

34 FORMAT (//lOX, 'ELEMENT CONNECT IVITY INFO. FOR ZONEZ'", 

!//, 6X, •" lELEM", 16X, "LOCAL NODE NOS",/, 

! 1 6X , ■" 1 ■" , 5 X , ■" 2 ' , 5X , 3 ' , 5 X , 4 , 5X , 5 ■ , 4 X , ' 6 ' , 6X , ' 7 ■ , 5 X , ' 8 ’ , / ) 
DO 510 IELEM-=1 ,NELeM2 

WRITE (24, 11) lELEM, (LN0DS2( lELEM, JJ) , JJ=1 , 8) 

510 CONTINUE 

DO 100 IELEM=1,NELEM1 


3 


n n o no 


I 


IM.l 1.50 IN(]DP =1 , NMUDP 

LMODO' JELEM, INODP ) -'LNDDSl <. lE-'L.EM, iriOOP 1 
150 COMTTMUE 
100 cowriMUE ■ 

m l p’ 1 ■ ^ MP 0 1 M 1 (. 2 wM X 1 + 1 ) 

NS I F2=NS I F 1 +2w ( NX 1 --MELCOM ) 

DO 220 IELEM-=1,NELEM2 
DG 250 INODP=l,NNaDP 
NBHI.FT>'NSIF2 

IF ( (. lELEM. LE. NELCOM) . AND. (. INODP . LE. 3) NSHir T-t'1SIFl 
LNODS ( I ELEM+NELEN 1 ^ I NODP .1 =LN0DS2 ( lELEM , INODP ) +NSH IFT 
Il-IELEM+NELEMl 
250 CONTINUE 
220 CONTINUE 

WR I TE ( 24 , 25 ) NEI.EM , NPO I N 
25 FORMAT (//5X, ■'NELEM=-', 15, 2X, • NPOIN- - , 15/) 

WRITE f. 24, 35) 

35 FORMAT (//lOX, 'ELEMENT CONNECTIVITY INFO. FOR DOMAIN', 

! //,6X, 'lELEM', 16X, 'LOCAL NODE NOS',/, 

! 16X, 'l ',fcX, •'2',6X, "3",6X, ■'4',6X, 'S^BX, •'6',7X, ■'7',6X, 'S',/) 
DO 520 IELEM=1 ,NELEM 

WR I TE (24,11) I ELEM , ( LNODS ( I ELEN , J J ) , J J= 1 , 8 ) 

520 CONTINUE 
RETURN 
END 


SUBROUTINE ZONE (NX , NY, NPQINX , NELEMX , LNODBX ) 
DIMENSION LWODSX( 111,8) 

DO 10 I =--^1, NELEMX 
DO 10 J^-1,8 
10 LNODSXd, J)”0 

DO 100 IELEM=1 , NELEMX 
IY="IELEM/NX 
Y--=FLOAT( lELEM) /NX 
D=^=Y I Y 

ir (D.PO.O.O) IY=dY“l 
IX-IELEM-IYwNX 

t.NODSX ( IEI,..,EN, 1 ) -IY>= (3hMX+2) + ( I X- 1 ) ^<2+1 
LNODSX ( I ELEM , 2) -LNODSX ( I ELEM , 1 ) + 1 
LNODS X ( I ELEM , 3 ) =^LNODSX ( I ELEM , 2 ) + 1 
LNODSX ( I ELEM , 4 ) ^---LNODSX ( I ELEM , 3) '-S^NX ~ I X M 
L..N(;!1)SX ( I EI..EM , 5 ) ■■■4.N0DBX ( I ELEM , 4 ) +NX + 1 X + 1 
' LNODSX ( I ELEM , 6 ) '-^LNODSX ( I ELEM , 5 ) - 1 
LNODBX ( lELEM, 7) ^LNODSX ( lELEM, 6) -1 
LNODSX ( lELEM, S)==LN0DSX ( lELEM, 4) -1 
100 CONTINUE 
RETURN 
END 


SUBROUTINE AUTO(NXE'L, NYEL, NBCON, XNODX , YNDDX , MTOT XY , I CALL) 
DIMENSION TRSP0L(2000) 

INTEGER BNfJD(lOOO) 
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D I rir.'.NS I (IN X NOD X (. 2000 ) , YrJOD X 2000 ) 

COriMOM/ AF>>EA2/NX , NY , NOD , NEL, NDXDIR, NDYDIR 

NX~NXEL 

MY-^MYEL 

MDXDIR^-2mNX (-1 

MI)YDI,Rr.-2icNY+l 

MBf:0N--4w(MX (-NY) 

ME.L="Y'1X^^MY 

(. ?-''-MX’( 1 ) ?>; ( 7^NY+1 ) - NEL 
CALL INPUT(NBCON,NOD, ICALL) 

CALL. BDYMOD (. NBCON , BNOD ) 

CALL INI POL (NDXDIR, NDYDIR, MBCOM,NTaTXY, BNOD, TRBPOL, ICAL.L.1 
DO 133 1=1, NBCON 
133 CONTINUE 
LAKE=1 

DO 53 IBM=1,N0D 

XNODX ( IBM.) = T RSPOLTLAKE) 

YMODX (IBM) -=TRSPDL(LAKE+1 ) 

LAKE=LAKE+2 
53 CONTINUE 

DO 1239 1=1, NOD 
1289 CONTINUE 
RETURN 
END 

C 

c 

SUBROUTINE INPUT (NBCONX , MODX , ICALL) 

INTEGER BNODX(IOOO) 

COMMON / A2 / BOUDC 1 (48,2), B0UDC2 (00,2) 

COMMON/ ARE A2/ NX , NY, NOD, NEL, NDXDIR, NDYDIR 
COMMON/ A4/ XO, YO, XI, Y1,X2,Y2 

C0MM0N/A5/RAKEAN, DCUT, CHPTHK , CONLEN, C03A, SINA 

NN0DP=--8 

13=7 

IF (ICALL. EQ. 1 ) GOTO 500 
IF(ICALL.EQ.2) GOTO 600 
500 CONTINUE 
11=0 
1 2"“ 29 
DX1=2. /6. 

DX2"^ DXl 

C' FIRST POUND C(3fJRD.S< SECOND 
XC0Rl=-3. 0 
YC0Rl=-0.6 
XC0R2=XC0R1 
YC0R2=0.0 
DO 10 I ="1,7 
I 1 1 1 'M 
12= 12-11 

BOUDC 1(11,1) =' XCORl +FLOAT ( I -1 ) wDX I 
BOUDC 1 (II, 2)=YCaRl 
B0(.)DC1 ( 12, 1. )="XCGR2-+FL0AT( I -1 )«DX2 
BOUDC 1 (12, 2)=YC0R2 
10 CONTINUE 
14 =36 



DX3 -2.0/12 
dX3 

XCf.)(’.3^-- -1.0 
YCC)n:3=--0.6 
XC(:)l’.4-^^XCOR3 
Yf.;0R4=-0. 0 
DO 20 
I3^---.T3-<- 1 
, I4---=T4»1 

BOUDC: 1(13,1 ^ =^XC0R3+FL0AT ( 1 1 fsDX3 
BOUDCl ( 13, 2.)--=YC0R3 
ED'JDC 1(14,1 ) =Xr.0R4 +FL0AT ( I ) sDX4 
B0UDC1(I4,2)=YC0R4 
20 CONTINUE 
J5==18 
DY5”^0. 125 
XC0R5L“-3.0 
YCOROL " -0 . 6 
XC0R3R-- 1 . 0 
YC0R5R=^YC0R5L 
DO 30 I =1,4 

BOUDC 1 ( J5 , 1 ) =XC0R5L 
BOUDC 1 ( J5 , 2 ) =YCaR5L+FL0AT ( I ) ><DY5 
, BOUDCl ( J5+1 , 1 )=XC0R5R 

BOUDC 1 ( JS+ 1,2) =YC0R5R+FL0AT ( I ) ?^DY5 
30 CONTINUE 
J6=-26 
DYA’^0.05 
XCUR6L=-3.0 
YC0RAL=“0. 1 
XCOPY-TU 1 . 0 
YC0R6R=-0. 1 
DO 40 1=1,1 
JA=J6+2 

BOUDCl (J6, 1)=XC0R6L 
BOUDCl (J6, 2)=YCaR6L+FLaAT(I)sDY6 
BOUDC 1 ( 3 6-+ 1 , 1 ) =XC0R6R 
BOUDC 1 ( J6-<- 1,2) =YC0R6R+FL0A Id) sDY6 
40 CONTINUE 

DO 37 KK=1,NBC0NX 
I7=KK 

37 CONTINUE 
GO TO 700 
600 CONTINUE 
1 1=0 
1 2=6 

DXl-=2.0/6.0 
DX2=1. 0/6.0 
XCORl =-3.0 
YCOR 1=0.0 
X COR 2 =-1.0 
YC0R2=-0.0 
DO 110 1=1,7 
11 = 11+1 
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1 T 2 M, 

Pri'JDC;? U 1 , 1 ) = XnORl+FLOAT ( I -1 ) >:1)X 1 
BOi..'OC2(Il,2)=YCORl 
EiDUDt'2 (12,1 )=^XC0R2+FLQAT ( I -1 ) 
B(3UDC2(I2,2.)~YC0R2 
110 CONTINUE 
1 3" 13 
NELRF=5 

NN0DRF=NELRFs2 

dl3^t:dmlen/nnddrf 

DO 120 I"1,NN0DRF 
13=1341 

B0UDC2 (13,1) =FLOAT ( I ) sDL3«9 I NA 
B0UDC2 (13,2) =FLOAT ( I ) wDL3wC03A 
120 CONTINUE 
14=23 

DL4='C0NLEN/6.0 
DO 130 I ^=1,6 
14=14+1 

B0UDC2 ( 1 4 , 1 ) = ( CONLEN+FLOAT ( I ) s^DL4 ) I N A 
B0UDC2 ( I 4 , 2 ) = (CONLEN+FLOAT ( I ) sDL4 ) sCOSA 
130 CONTINUE 

X 2=2 . OwCONLENwS I NA 

Y2= '2 . 0-^C0NL.EN>*C0SA 

X1=X2 ■CHRTHKmCOSA 

Y1=Y24CHPTHKwSINA 

J51=28 

J52=40 

DUr'il “O.O/A.O 
DL52- 0.50/6.0 
X COR 5- -3.0 
YC0R51="0.0 
YC0R52=0 . 5 
J61=29 
J62=41 

' I)L61=0. 3«CHPTUK/6.0 

DL62 =0. 7wCHPTHK/6. 0 
XCORAI-X? 

YC0R6 ] "-Y2 

XC0R62=X2 -6 . 0«DL61 wCOSA 
YC0R62=Y2+6 . 0»DL6 1 wS INA 
DO 140 1=1,6 
J51=J51+2 

FinUDC2 ( J51 , 1 )=XC0R5 

B0UDC2 ( J5 1 , 2 ) =YC0R5 1 +FLOAT ( I ) DU5 1 

J52="J524 2 

B0UDC2 ( J52 , 1 ) -'XCORS 

BOUDC? ( J52 , 2 ) =YC0R52.+Fl..0AT ( I ) «DL6 1 ^sCOSA 

vJ6C“J61+2 ' 

P0UDC2 ( J61 , 1 ) =XC0R61 -FLOAT ( I ) «DL61«C03A 
BaUDC2( J61 , 2 )=YC 0 R 61 +FL 0 AT ( I ) sDL6 1 -.^SI NA 
I J62=J62+2 

BOUDC? ( J62 , 1 ) =XC0R62 "FLOAT ( I ) S'DL63'KC0SA 
B0UDC2 ( J62 , 2 ) =YC0R62+FL0AT ( I ) sDL62?^S I NA 
140 CONTINUE 
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CALI, CURVF ( PI.'CDC? '> 
' DO 36 , NBf'.ONX 

I 7--KK 

36 CrjHTINUE 
700 RETURN 
EAJD 


SUBRO'JT I NE I MTPni . ( NNX , NNV , T!RCON , r.!T OT X Y , DMOD , TPRPni., , I CALL ) 
DTMFM3I0N TR3P0L (2000) , BOUDV(800, 2.) 

INTEGER BNOD(IOOO) 

COMMON/ A2/B0UDC1 (48,2), 800002(80,2) 

COMMON /AREA2 /MX , NY, NOD, NEL. NDXDIR, NDYDIR 
D T MENS I ON X A 1 ^ 500 , 500 ) , X A2 (500,500 ) 

D I MENS I ON YAK 500 , 500 ) , YA2 ( 500 , 500 ) 

IF( ICALL.EO. 1) GOTO 30 
IF(ICALL.En.2> GOTO 40 
! 30 CONTINUE 

DO 31 I=1,NBC0N 
BOUDV (1,1) ■••=B0UDC 1(1,1) 

B0UDV(I,2)=B0UDC1(I,2) 

31 CONTINUE 
GO TO 50 
40 CONTINUE 

DO 4 1 I ), , MBT'.ON 
BOUDV (1,1) ^■•-B0UDC2 (1,1) 

BOUD V ( 1 , 2 ) =^B<0UDC2 ( 1 , 2 ) 

' 41 CONTINUE 
SO CONTINUE 

DO 450 I = 1,NBC(3N 
DO 450 J=-l,2 
450 CONTINUE 

DO 1551 I01=--1,NDXDIR 
K44 = l 

K55^^2^;NX+K44 

DO 1552 JOl'^KNDYDIR 

RI 1=FL.0AT ( 101 -1 ) /FLaAT.(NDXDIR-l ) 

RI2=FL0AT(NDXDIR-I01)/FL0AT(NDXDIR-1) 

IF( JOl -EO.NDYDIR) KSS^NBCON 

X A 1 ( 1 0 1 , JO 1 ) =R 1 1 mBOUDV ( K55 , 1 ) +-R 1 2?^B0UDV ( K 4 4 , 1 ) 

Y A 1 ( 1 0 1 , JO 1 ) =R 1 1 wBOUD V ( K55 , 2 ) +R 1 2^4 BOUDV ( K 44 , 2 ) 

K44=K55+1 
K55=K44+1 
1552 CONTINUE 
1551 CONTINUE 
K99"l 

K77”NDXDIR i'2?< (NDYDIR-2) + 1 
DO 1777 I=1,NDXDIR 
DO 1778 J=1,NDYDIR 
R J 1 ;=FLOAT ( J - 1 ) /FLOAT ( NDYD I R - 1 ) 

RJ2=‘T- LOAT (NDYDIR-J) /FLOAT (NDYDIR--1 ) 

X A2 ( I , J ) .^R J 1 w ( BOUDV ( K77 , 1 ) -X A 1 ( I , NDYD I R ) ) +R J2w ( BOUDV ( K99 , 
! 1)-XA1(I, 1)) 

T I P-^-BOUDV ( K77 , 2 ) -YA 1(1, NDYD I R) 

TEK- BOUDV (K99, 2) -YAl (1,1) 
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YA2 ( I , J) -^^RJ 1 wTIP fRJ2»TEK 
1778 CONTI MUR 
K 77^0 <77+1 
K99r.K99+l 
1777 CONTINUE 
MM2- 1. 

DO 2300 J^-1 ,NDYDIR,2 

DO 2330 I^1,NDXDIR 

TRSPOL <.MM2V--= XA1 < I , J) +XA2 < I , J) 

LMN2“MM2+1 

TRSPOL LMM2 ) =YA 1 ( I , J ) +YA2 ( I , J ) 

MM2^=^MM2+2 
2330 CONT I NUE 

MM2^^■=MM2■^2^o^MX+2 
2300 CONTINUE 

LTTL=-Y2w<.NDXDIR) +l 

DO 2312 NDYDIR-1,2 

DO 2310 I "1 , MDXDIR, 2 

TRSP OL f. I TT 1 ) =X A 1 (. I , J ) +XA2 < I , J ) 

MLrr-LTTl + 1 

TRSPOL (. ML T. T I =:YA 1 ( I , J .) +YA2 ( I , J ) 
LTTl“LTTl+2 
231'0 CONTINUE 

LTTO=l.,TTl+2»NDXDIR 
2312 CONTINUE 
RETURN 
END 

C 

c 

BUEROUT I NF. BDYNOD (. NBCQN , BNOD ) 

INTEGER BNOD (.1000.) 

C0MM0N/AREA2/NX, NY, NOD, NEL, NDXDIR, NDYDIR 
LIV"0 

DO 1880 J4=^l, NDXDIR 

LIV™LIV+1 

BNODd .IV)==J4 

LACK =--ND X D I R+2w ( NDYD I R -2 ,) + 1 
MA55=N..IV+LACK~1 
I BNaD(MA55.)=N0D-NDXDIR+J4 
1880 CONTINUE 
IC0UNT::=1 
KE.E:T“H NDXDIR 
DO 2777 JAT=2,NDYDIR-1 
BNOD ( K EET ,) “BNOD ( K EET - 1 .) + 1 
IAEA ■KEET+l 

BNOD ( I ASA,) ='BNOD (KEET ) +NXsICOUNT 

ICOUNr==ICOUNT+l 

KEET “KEET +2 

■ IF( icourir.EQ.s.) icount^^i 
2777 CONTINUE 
RETURN 
END 

C 

SUBROUT’ I ME CURVE ( B0UDC2 .) 
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D I nrr I'" I ( ir-j f3gui')(J 2 ( ao , 2 ) , anglm ( aaa ) 

CC.U'ini 1!'l / A'1 /VO, YO , X 1 , Y 1 , , V-^ ^ Y 7 

C0MM0N/A5/RAKEArj, DCUT, CHPTMK, CCMLEN, CORA, RTMA 
DO 200 ,380 

ANGLOK T .''=0. 0 

200 continue: 

SLOP^-COSA/SINA 
PI=- 3. 1 ^ 159.77 
X0--3.0 
Y0=- - 1 . 05 

B0UDn3C52r 1 .1 ^XO 
E''n'.'OC? U.';7, 2.''--=Y0 
AI..FA^-(Yt "I. 01 / (YO -1 .0,1 

YBAFT ( ( YO - AI..F A>-Y 1 ) - BLDP^^ (. XO -ALFA^^ X 1 ) ) / ( 1 . 0 -ALFA ) 
C ( Y 1 "■GLOP wX 1 -YBAR .) w (. Y 1 - 1 . 0 ) 

XL.EhL^^-X1. ■'•XO 

DX-XLEM/IO.O 

DXl-"3-0.?vDX/A.0' 

DX2 ■3)X1 

DX3-"3.0;<^DX/10.0 
DX'LOIX/A. 0 
XOL-XO 

X03-^^X01 f-3. 0'>-DX 
X03' X02'i 3. 0.".:I)X 
X0^-^'^"X03'-3.0«DX 
IB0UDL'--52 
I BGUD2"Y'5a 
IB0UD3''-A^ 

IB0UD'V--7T 

' Il^"'35Ei! I 

J1’Y36'T 
Kl=''370 
L 1 ^^Y380 

C EON : (. Y -l-IX -■ YBAR ) >;■ ( Y - 1.0) == ( Y 1 -MX 1 -YBAR ) (. Y 1 - 1.0) 

C SLOPE FROM (FOM: OYDX^^M ( Y- 1 . 0) / ( (, Y-1 . 0) + ( Y -MX -YBAR ) ) 
DO 100 I ^--1,10 
IF( I.Gr.6) GOTO 23 
XCORL- XOH InOXl 
RL^-^SLOPir-XCORl fYBAR 
Cl='ni ■0 

YCCtRl--(.Rl“i LinORK (.Rl + 1 ) «-s2-4 . 0?^C1 ) )/2.0 
D Y D X Sl... OP w ( YCOR 1 - 1 . 0 ) 

DYDX="T)YDX/ ( ( YCOFII -1 . 0) + (YCORl -SLQPwXCORl -YBAR) ) 
BETA'^^ATANIDYDX) 

ANBI, E-^ 90 . 0 - BFTA-wl 80. 0/PI 

I .rBOUD^■^TBaUDl + l 
I1==I1"H 

BOUDCLK IPOUDl , 1 )=-^XCGRl 
B0UDC2 (. I BOUD 1,2) =YCOR 1 
ANBLM ( 1 1 ) ==ANGLE 
XC0R2--X(:)2'MwDX2 
R2-3l„GR XCGR2+YBAR 
CL-=R2-C 

YCC:)R2=- ( R2-'' 1 + SORT (. i R2+ 1 ) -4 . 0«C 1 ) ) / 2 . 0 

DYDX-=3L0P«-<,YC0R2-1 .0) 


10 


rvVnX -OYDX/ ( (YCGR 2-1 . O) +-(Yf:GR?. -si .OP^r-XCURlZ ■YBAP,.'' ) 

BE i A-ATAN(DYUy,) 

AMGl R ^^^70 . 0 -eE TAm 1 80 . 0/P I 

J,BnUD2=IB0UD24'l 

J1:'=J1 + 1 

E^DUDtX.':-! ( I BL)UD2 , 1 ) =XC0F12 
B0UDC2 ( I BQUD2 , 2.^ =YC0R2 
AMGLN( J1 ) WANGLE 
>:C0RA'=X04-<-TwDX4 
R '!■ S !. . 0 P « X C [3 R4 +• Y BAR 
CB-RA-C 

YCOR 4 i R4-I- 1 +SnRT (. (. R4+ 1 ) ws2 -4 . OsC 1 ) ) / 2 . 0 
DYDX' R! ,.0P^;(YC0R4”1.0) 

DYDX ^ DYDX/ F rt'CQR4-l . 0) -i- ( YCDR4-5L0P^tXC0R4 -YBAR) 3 
BR"rA-^ATAN(DYDX) 

AMBL r-:='90. 0-PETA^<<180. 0/PI 
IBOUD4'-IBGUD4 + l 
L 1 =■ i . .1. -I- 1 

BUUBC7 t IBGUD4, 1 )=XC0R4 
BOUnC? (. I BGUD4 , 2 3 ==YC0R4 
AMGLM'.Ll 3-^ ANGI„E 
23 CGNTTNUE 

XCGRR •■.X034 I;f4:)X3 
R3-SL0PWXC0R3 i-YBAR 
C1,“=R3'4:; 

YC'GRR" ( R3 (• 1 -I RGRT ( ( R3+ 1 3 «s2 -4 . O^C 1 3 3 / 2 . 0 
D YD X ■•■<5L0P ( YCGRB -1.03 

DYDX-DYDX/ ( ( YCC.)R3- 1 . 03 + ( YCC)R3-SL0P^^XC0R3 “YBAR3 3 
' BETA^ArAN(DYDX3 ' 

ANGLE=-^^90.0-EETAwl80.0/PI 
IB0UD3=IB0UD3+1 
K 1 =■ 1 + 1 

B(::iUl)C2 ( I BGUD3 . 1 3 =XC0R3 

B0UDC2(IB0UD3,23=YC0R3 , 

AMGLNCKl 3=ANGLE 
100 CONTINUE 

DO 234 ,386 ' 

WR I TE f 27 r w 3 I , ANGLN C 1 3 
234 CONTINUE 
RETURN 

END . i 


1 1 




n a n n n 


c FINITE ELEMENT ANALYSIS OF 

C tttttttttttttttt+IttttttttWtttt ORTHOGONAL METAL CUTTING tttttttttttt 

C 

C MASTER FLUID 
C 

IMPLICIT RFAL^^a(A-H,0-Z) 

COMMON/ARCl /POSGP (3) , WEIGP (3) 

C0MM0N/ARC2/LN0DS( 1 11,8), COORD (386, 2) . NADFM(386) , N0DFM(386) 
I COMMON/ ARC3/VARB1 (2000) , VARB2( 2000) 

COMMON/ ARC35/EPSN ( 1 1 1 , 9 ) , XG ( 1 1 1 , 9 ) , YG ( 1 1 1 , 9 ) , V I SCP (111,9) 

! , DAREAG (111,9) 

COMMON/ ABC/ I R 

COMMON / ARC99 3 / T REF , P ECLET , CONR AT , T FBT AR , HST AR , 
!BSTAR,SIGYO,SPHTO,THKO, SIGMA 
DIMENSION EQRHS(2000) 

D I MENS 1 ON LPOUD ( 2000 ) , LHEDV ( 1 19), PNORM (119), EOUDV ( 2000 ) 

D I MENS I ON GFLUM (119,119) 

C 

CHARACTER «00 FNAME 

WRITE(w, w) -ENTER INPUT FILE NAME" 

READ ( w , •■■ ( A50 ) ■' ) FNAME 
OPENIUNIT - 20, FILE"FNAME) 

OPEN ( UN I T ==40 , F I LE= •' VAR . I NP ■- ) 

OPEN ( UN I T=-=24 , F I LE= ■' S3V3TS . OUT - ) 

0PEN(UMIT^35, FILE= 'S3V3VTS. OUT- ) 

OPENIUNI T~36, FrLE=-S3V3TTS. OUT- ) 

0PE:N ( UN I T-28 , F I l.E- - S3V3ETE! . OUT ■ ) 

OPEN ( UN I T=88 , F I LE^^ •- TEMPTS . DA T - ) 

OPEN ( UN I T>=90 , F I LE= ' STRA I NTS . DAT •' ) 

OPEN ( UN I T=38 , F I LE= ■' STRESTS . DAT •' ) 

C 

SET DYN. DIM. VALS 
READ (20, w) IR 

CALL DIMENS(MELEM,MFRDN,MPOIN,MTOTV) 

READ ALL PROB. DATA 
CALL 1) I NPUT ( I AXBY , BOUDV , LBOUD , MELEM , 

! MPOIN, MTOTV, NBCON, NDOFM, NELEM, NEVAB, NGAUS, NITER, 

! NNODL, NNODP , NPOIN, NTOTV, RELAX, TOLER, 

! XFORC, YFORC) 

1176 CAI.L ITERAT ( I AXBY, GFLUM, BOUDV, LBOUD, LHEDV, PNORM, EDRHS, 

! MELEM, MFRON, MPOIN, MTOTV, NBCON, NDOFM, 

! NELEM , NEVAB , NGAUS , N I TER , NNODL , NNODP , NPO I N , NT OTV , 

! RELAX , TOLER, XFORC, YFORC) 

WRITE (28, 100) 

WRrTE(w, 100) 

100 FORMAT (20X, --STRAIN RATE DISTTN'-//, 

! -'KELEM-TASX, 'GAUSS POINTS-',/, 

! 18X, 1 - , 12X, '2' , 12X, -'3-', 12X, --4-', lOX, 'S', lOX, ■-6', lOX, -'7' , 

' 12X -'8'- 12X --9' /) 

100 FORMAT (// 'EL E:M G.PT X -COORD Y- COORD STRAINRATE ' , 

DO 200 1. MELEM 

DO .1 1 , 9 


I 


I 


1 


WRI TE (28, 1^.) I , J, XB( I , J) , YGd , J) , EPSM( I , J) 

C 111 F0RMAT(I3,2X, I2,D8.4,2X,D8.^,2X.D10.4) 

ENDDO 

200 CONTINUE 

CALL FORCE 
CALL STRATES 
CALI- CHTHK 
TYPE'w, -MAIN OVER- 
STOP 
END 

^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ kf; ^ 

c 

BUPROUT I NE D I MENS ( MELEM , MFRON , MP 0 1 N , MTOT V ) 

IMPLICIT REAL^^aiA-HrO-Z) 

ME.l,.ETL = 1 1 1 
MFRGN-d 19 
MP0TN"306 
M TO 10^^2000 
RETURN 
END 

C/ ^ i'f: ^ 'S^< -p- ^ ^ ^ ^ ^ ^ ;Ct ^ ^ ^ ^ X. ^ ^ ^ ^ 

SUBROUTINE DINPUT( I AXSY, BOUDV, LBOUD, ' 

! MFLEM, MPDIN, MTOTV, NBCON, NDOFM, NELEM, NEVAB, 

! MGAUS , N I TER , NNODL , NNODP , NPO I N, NTOTV , RELAX , 

I TOL ER , XFORC , YFDRC ) 

IMPLICIT REALw8(A-H,0-Z) 

C 

C0MM0N/ARC2/LNDDS(111,8) ,CnaRD(386,2) ,NADFM(3a&) ,NDDFM(386) 
DIMENSION LBOUD ( 2000 ), BOUDV (2000) 

C0MM0N/ARC3/VARB1 (.2000) , VARB2(2000) 

COMMON/ ARC 1 0/CUTVEL , DCUT , W I D , RAKE , GAMMA 
COMMON/ARC 12/FRIFAC,P I, EN,C0SA,SINA,SIN2A 
COMMON/ ARC 1 3/ RHO,MBOUEL 

CQMM0N/ARC14/IBEL(.60) , I SURF (60) , I BIDE (60) 

COMMON / ARC30 / ANGLN ( 386 ) 

COMMON /ARC991 /TREE, PECLET, CONR AT, TFSTAR, HSTAR, 

! BSTAR , S I GYO , SPHTO , THK 0 , S I GMA 
' COMMON /ABC/ I R 

C 

RNO^'/OBO.O 
NDOFM "4 
MEVAB=^^2S 
NGAUS-^3 
MM0DL..=4 
NNODP ^8 

READ (20, 1000) TITLE 
1000 FORMAT! 12A6) 

READ(20,^<*) lAXSY, NELEM, NITER, NPOIN,NRPDN 
READ ( 20 , M ) N I CON , NBCON , NEBCN , NNBCN 

CALL. DI AGMl (NBCON, NEBCN, MELEM, NICON, NNBCN, MPDIN, NRPON) 

READ(20, w) RELAX, TOLER, XFORC, YFORC 

READ ( 20 , ) EN , GAMMA , B IGYO 

READ ( 20 , w ) CUT VEL , DCUT , RAK E , P'R I FAC , W I D 

PI ^’"3. 1415927 

RAKFAM -TIAKEwP I / 180 . 0 


RAKE2^2.0;^RAKEAN 
’ 9IN2A“DBINf.RAKE2) 

COS A--DCOS ( R AK E AN ) 

BINA=^'1.)SIM(RAKEAN) 

P'TI AFR- (CUTVEI ./ ( 1, . 732-<t:DCUTwRAMMA) ) ( 1 . O/FFO 

BSrAR-BSTAR/GIRYO 

READ (.20 ) FI , T F1K0 , THKTOL , RHO . SPHTO 

READ(20,w) TAMB, TWP 
AF. P FI A^^'TFIK 0 / ( FU in=»;BPNTO ) 

P ECI. J-. T ^=CU T VE I DCUT / ALP HA 
TRFF=^BIGYO/ (RHOwBPFnO) 

HSTAR=HwDCUT/THKO 
TFSTAR"TAMB/TREF 
CONRA T^^TI IKO/ ( THKTQL+THKO) 

BIRMA^^BIf^YO/ (RH0wCUTVELw«2) 

WR I TE ( 30 , 355 > DCUT , CUT VEL , R AKE , H , THKO , SPHTO , RHO , FR I FAC 
355 FORMAT f.//3X, DATA USED', 

! /2X, 

! //4X, -'DEPTH OF CUT F12. 8, 3X, ■CUTTING VEL.: ■",1)1 

! /4X, "RAKE ANGLE : ' , D12. 4, 3X , "HT. TR. COEFF. ; ' , D1 

! /4X, "THERM. COND. : , 1)12. 4, 3X, 'SP. HEAT .-^Dl 

! /4X, -"DENSITY : " , D12. 4, 3X, "FRIC. FAC. ;',D1 

! /3X, ) 

DO 10 IPOIN^l ,NPOIN 
DO 10 I DIME" 1,2 
COORD! IPOIN, IDIME.)=0.0 
10 CONTINUE 

DO 20 IPaiN-l,NRPON 
JPOIhF-IPOIN 

RFAD(20,-<«) JPOIN, (COORD! IPDIN, IDIME) , IDIME=1 , 2) 

20 CONTINUE 

READ ! 20, SCALE 
DO 25 IDIME=1,2 
DO 25 IP0IN=^=1,NRP0N 

COGR D ! I P 0 1 N , I D I ME ) =C00RD ! I P 0 1 N , 1 1) I ME ) »SC ALE / DCUT 
25 CONTINUE 

DO 30 IELEM"M,NELEM 
JELEM-=IELEM 

FIE AD!20, *4) JELEM, (LNODS! JELEM, INODP) , INODP-l , NNODP ) 

30 CONF'INUE 

READ ! 20, ^f.) MBURF 

READ !2C), ^<^) MFLBl , MELB2, MELB3, MELB4, MELB5, MELS6, MELS7, MEL58 

MB0UEL~MEL31 •<-MELS2+MELS3+MELS4 ^NELS5■^-MELS6+MELS7+MELS8 

DO 9 II^^1,MPDUEL 
J J " ” I T 

READ ! 20, .) J J , I SURF ! 1 1 ) , I BEL ! 1 1 ) , I BI DE ! 1 1 ) 

, 9 CONTINUE 

SET UP A VECTOR GIVING D.O.F AT EACH NODE 
ITEMP'=NN0DP-1 
DO 40 IELEM=1,NELEM 
DO 40 IN0DP==1, rrEMP,2 
NODFM ! LNODS ! I ELEM , I NODP ) ) "MDOFM 
NODFM ! F.NODS ! I ELEM , I NODP 1 ) ) =-=NDOFM - 1 
40 CONTINUE 


w w w w 


no no n o on 


c INTERPOl.E MIDSIDE NODE NHEP.E NOT KNOWN 
GO TO 71 

DO 70 lELEM^^l , NELEM 
DO 60 IN0Dr'-:2,NN0DF\2 
I PO I N =LMGDS ( I ELEM , I NODE ) 

1 EMP Y” DABS ( COORD (I PO I N , 1 ) ) +DABB ( COORD U PO I M , 2 1 ) 

I F ( TEMPY . NE . 0 . 0 ) GOTO 60 
JPO I N^L,,NODS (. I ELEM , I NODE - 1 ) 

KN0DP=^IN0DP + 1 
IF (KNODP . GT . NNODP ) KN0DP=1 
KPO I N'^LNODS ( I ELEM , KNODP ) 

DO 50 IDIME^1,2 

COORD ( I P 0 1 N , ,I D I ME ) = ( COORD C JP 0 1 N , I D I ME ) +CDORD ( K PO I N , I D I ME ) ) »0 . 
’50 CONTINUE ^ 

60 CONTINUE 

70 CONTINUE 

71 CONTINUE 

SETUP A VECTOR WITH FIRST D.O.F AT EACH NODE 
NADFM U ) = 1 
DO 80 IPniN^2,NP0IN 

NADFM ( I P 0 1 N ) “NADFM ( I PO I N - 1 1 +NODFM ( I PO I N - 1 1 
80 CONTINUE 

NTOTV:" NADF M (. NP 0 1 N ) +NODFM f. NP 0 1 N ) - 1 
INITIALISE REMAINING ARRAYS 
DO 90 IT0TV==1,|\1T0TV 
LBOUDC rT0TV)=0 
EDUI)V(ITOTV):^0.0 
VARBHrr0TV.')=0’.0 
I VARB2(IT0TV)=0.0 
90 CONTINUE 

IFn,R.EO. DTHEN 
DO JT0TV=1,NT0TV 

READ<.40,«)IP0IN, IDOFM , VARB2 ( JTOTV) 

END DO 

DO UC0N--1 ,NICON 
RE AD (20,=^ I 
ENDDO 

GOTO no 

END IF 

wwwwREAD INITIAL CODN.S IF ANY 
IF(NICON. EQ.O) GO TO 110 
DO 100 IICaN=n,NICON 
READ ( 20 r-*) IPOIN, IDOFM, TEMPY 
IF ( I DOFM. GT . 1 ) I DOFM=NDDFM ( IPOIN.) 

JTOTV“NADFM ( I PO I N ) + I DOFM - 1 
VARBl (.JTOTV)=TEMPY 
VARB2 ( JTO TV .) =TEMP Y 
100 CONTINUE 
110 CONTINUE 

CHECK NODAL CONN.S 8. COORDINATESX 

CAIJ.. D I AGN2 ( MELEM , MPO I N , NELEM , N I CON , 


c 


!NNnDP,NPOIN,NTOTV) 


C 

C 


C 

c 


DO 1 70 I BCON^' 1 , NBCOM 

READ(20,w) IPOIN, IDaFM,BVAIJJ 
IF((IPDIN.EO. D.AND. (ID0FM.EQ.2))-fHEN 
WRITE (35, 799) BVALU 

799 F0FWAT(/5X, -PRESSURE AT NODE NO. 1 IS •',08.2, ' OF SIOMAY') 
EWDIF 

I F ( I DOFM , EO . 4 ) EVALU=BVALU/T REF 
rTOTV==NADFM( IPOIN) + ID0FN-1 
IF(N0DFM(IP0IN).EQ.4) GOTO 238 
I F ( ( I DOFM+ 1 ) . GT . NODFM ( IPO IN) ) I TOTV- 1 TOTV -1 
238 LBOUD( IT0TV)=1 

BOUDV ( I TOTV) =BVALU 
60 TO ( 1 40,1 50 ,160,165) I DOFM 
140 CONTINUE 
GOTO 170 
150 CONTINUE 

GOTO 170 

IF(NODFM(IPOIN) .NE. (IDOFM+1) ) GOTO 170 
STOP 

160 CONTINUE 
GOTO 170 
165 CONTINUE 

JTOTV^ITOTV 
170 CONTINUE 
175 CONTINUE 

DO 990 IT0TV=1,NT0TV 
JTOTV"-.ITOTV 
990 CONTINUE 

DO 999 I=1,MP0IN 

READ ( 20 , « ) J , AN6LN ( I ) 

999 CONTINUE 

DO 400 KKL.= 1,MP0IN 
LL^--=KKL 
400 CONTINUE 
RETURN 
END 

ib; ^ ^ ^ ^ :fir. ^ *0" io; ^0; *0: ^ ^ tt): *0^ 


SUBROUTINE DI AGNl (NBCON, NEBCN, NELEM, NICON, NNBCN, NPOIN, NRPON) 
IMPLICIT REALw8(A-H,0-Z) 

DIMENSION NECUOdSO) ,NER0R(8) 

SCRUTINY OF CONTROL DATA 


DO 10 IFROR=^U. ,8 
10 NERORC IEROR.C-0 
SCRUTINISE CONTROL DATA S. 
IF(NPOIN.LE.O) 

IF (NRPON. LE. 0. OR. NRPON. 
IF(MFLEN,1..E. 0) 

I F ( NPO I N . GT . NELEMwB ) 

I F ( NBCON . GT . NP 0 1 Nw4 ) 
IF(NIC0N.GT.NP0INw4) 
IF(NEBCN.GT. NELEM) 


PRINT MESSAGE 

NER0R( 1)^=1 
GT.NP0IN)MER0R(2)=2 
NER0R(3)=3 
NER0R(4)==4 
NER0R(5)"5 
NER0R(6)=6 
NER0R(7)=7 


n n n n n n 


MERnR(a.v-e 


I F ( NMBCN . GT . NPO I M?=^4 .1 
jrFUI(>: 0 

DO 20 IER0n:=l,3 
I F (. MERDR ( I EROR ) , EQ . 0 ) GDI 0 20 
JER0R^=1 

WRITE (24, 2000) IFROR 
2000 FORMAT (/,33H CONTROL PARAMETER ERRORsw^sERROR, T5) 

20 CONTINUE 

IF(JEROR.EQ.O) RETURN 
WRITE (24, 20 10) 

2010 FORMAT (/ lOX, 37HDATA FOLLOWING ERROR IN CONTROL CARDS/) 
30 CONTINUE 

READ (20, 1000) NECHO 
1000 FORMAT! 150A1) 

WRITE (24, 2020) NECHO 
2020 FORMAT! 1 OX, 80A1) 

GOTO 30 

RETURN 

END 


' BUBROUT INE DIAGN2(NELEM, MPOIN,NELEM, NICON 
! ,NNODP,NPOIN,NTOTV) 

IMPLICIT REAL«8(A-H,0-Z) 

D I MENS I ON NECHO (150), NEROR (13) 

C0MM0N/ARC2/LN0DS(1 1 1,8), COORD (386, 2) , NADFM(386) ,N0DFM(386) 
DO 10 IER0R=1,13 
NEROR! I EROR) =^=0 
10 CONTINUE 
C CHECK PHY. PROPS 

I F ( DENS . !.„E .0.0. OR . V I BCY . LE .0,0 ) NEROR ( 9 ) == 1 ' 

C CHECK F'OR I DENT. COORDS 
DO 40 IPOIM'“2,NPOIN 
JPOITT^IPOIN-1 
DO 30 KP0IN==1 , JPOIN 

I F ( COORD ( I P 0 1 N , 1 ) . NE . COORD ( K P 0 1 N , 1 ) . DR . COORD ( I P 0 1 N , 2 ) . NE . 
!COORD(KPOIN, 2) ) GOTO 20 
NEROR! 10)- 1 
20 CONTINUE 
30 CONTINUE 
40 CONTINUE 

C CHECK ELEMENT NODAL NUMBERS 
C A -REAP I TAT I ON OF NOS. 

C B-NO. OUTSIDE PERM. BOUND 
DO 50 IP0IN=1 ,NP0IN 
DO 50 1ELEM=1,NELEM 
LC0NT=-0 

DO 50 IN0DE=^1,NNGDP 

IFd.NODS! lELEM, INODE) -NE. IPOIN) GOTO 50 
LCONT^LCONT< 1 
IF(LCONT.GT. 1 ) NEROR! 11 )-=l 
50 CONTINUE 

DO 60 lELEM^l ,NELEM 
DO 60 INODE"! ,NNODP 

I F ( LNODS ( I EL.E M , I NODE ) . LE . 0 . DR . LNODS ( I ELEM , I NODE ) . GT . NPOI N ) 


6 


n n 


iNERGRf. 12.)=1 
60 COM TIMUR 

C CHECK MO OF IMITIAL C0NDM3 

I F ( M I CON . GT . NTOT V ) NERDR ( 1 3 ^ 

JEROR^^^O 

DO 70 IFROR-IO, 1,3 
I F (. MFROR ( I EROR ) . EO . 0 ) GO FO 70 
JEROR^^ 1. 

WR I T E ( 2 A , 2000 ) I EROR 
2000 FORMAT(/10X, lOH DATA ERROR, 15) 

70 CONTINUE 

IF( JEROR.EQ.O) RETURN 
WRITE f.24, 2010) 

2010 F0RMAT(/10X, 14HREMAINING DATA/) 

80 CONTINUE 

READ (20, 1000)NECHO 
1000 FOFTMATdSOAl) 

WRITE (25, 2020) ECHO 
2020 FORMAT (/10X,80A1) 

RETURN 
END 

C 

SUBROUTINE ITERAT ( I AXSY, GFLUM, BGUDV, LBOUD, LHEDV, PNDRM, EORHS 
! , MFLEM , MFRON, MPO I N , MTOTV , NBCON, NDOFM, 

! NELEM , NEVAB , NGAUS , N I TER , NNODL , NNDDP , NPO I N , NTOT V , 

! RELAX , TOLER, XFORC, YFORC) 

IMPLICIT REALw0(.A-H,0-Z) 

C EXT. SUBROUTINES 
C PRESCR 
C FRONTS 
C WRITER 
C TOLREL 

COMMC IN / ARC 1 / P OBGP ( 3 ) , WE I GP ( 3 ) 

C0MM0M/ARC2/LN0DS( 1 1 1 , 8) , COORDS (336, 2) , NADFM (386) , N0DFM(.386) 

, C0MMGM/ARC3/VARB1 (2000) , VARB2. (2000) , 

C0MM0N/ARC35/EPSN( 1 1 1 , 9) , XG( 1 1 1 , 9) , YG( 1 1 1 , 9) , VISCP (111,9) 

! ,I.)AREAG(11 1,9) 

D I MENS I ON EORHS ( 2000 ) 

DIMENSION LBGUD(.2000) , LHEDV( 119) , PNORM( 119) , B0UDV(2000) 
DIMENSION GFLUM(119, 119) 

C SETUP ITER COUNTER ' 

IITER“0 

DO 1 I--=1,MELEM 
DO 1 J=l,9 
EPSNd , J)“0.0 
1 CONTINUE 
10 CONTINUE 

IITER“I I TER-t-l ’■ ' 

WRITE(35,w) •' IITER=--MITER 
DO 20 ITOTV-l ,NTOTV 

, e:qrhs(itotv)~o.o 

'20 CONTINUE 

C CALL FRONTS TO SETUP AND SOLVE GOVERN I G EONS, 

DO 50 Id, MFRON 


n o 


DO 50 J^^n,MFROM 
Bri JJM ( I , J ) ^ 0 . 0 
50 CONTINUE 
KOUNT^l 
NESPBC-0 

CALL FRONTS ( I AXSY, I ITER, MELEM, NFRON, MPOIN, MTDTV, NBCON, NELEM, 

! NEVAB, NNODL, NMODP , NPOIN, NTOTV, XFORC, YFORC, N0AUS, NESPBC, 6FLUM, 
! BOUDV , LBOUD , LHEDV, PNORM, EQRHS ) 

K0UNT^=K0UNT + 1 

CALL WRITER TO D/P ITER. RESULTS 

CALL WR I TER ( I I TER , NPO I N , NTOTV , NDOFM , NPO I N , HELEN) 
K0UNT^4<0UNT+1 

C CALL TOLREL TO CHECK CONVERGENCE AND RELAX VALUE IF NOT 

CAL.L TOLREL ( 1 1 TER, NTOTV, NCONV, NTOTV, RELAX , TOLER, VARBl , VARB2 
! , NADFM, MODFN, NPOIN) 

KDUNT===K0UNT+-1 

C RETURN TO N ASTER NO. IF ITERATIONS EXCEED NAXINUN 
IF (NCONV. EQ. 1) GOTO 12 
IF (IITER.LT. NITER) GOTO 10 
WRrTE(T«,2000) 

WRITE (35, 2000) 

2000 FORMAT (//32H SOLUTION HAS FAILED TO CONVERGE) 

12 CALL WRITERl IITER, MPOIN, NTOTV, NDOFN, NPOIN, HELEN) 

RETURN 

END 


C^5: ’.0: ^ iQ: ^ -c* >0: ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ 

SUBRDUT I ME WR I TER ( 1 1 TER , NPO I N , NTOTV , NDOFM , NPO I N , HELEN) 

I MPL I C I T REAL.?*8 ( A -H , 0 -Z ) 

DIMENSION VARBX(2000) 


COMMON / ARC2 / LNODS (111,8), COORD ( 306 , 2 ) , NADFM ( 386 ) , NODFM ( 386 ) 
C0MM0N/ARC3/VARB1 (2000) , VARB2(2000) 

COMMON/ ARC 1 0/CUTVEL , DCUT , RAKE • 

COMMON/ ABC/ I R 

C0MM0N/ARC99 1 /TREF, PECLET, CONRAT, TFSTAR, HSTAR, 

! BST AR , S I GYO , SPHTO , THKO , S I GMA 
WRITE (w, 111) I ITER 
111 FORMATdX, ••'IITER-^', 15) 

WRITE (35, 2000) I ITER 

2000 F0RMAT(//29H RESULTS FOR ITERATION NUMBER, 12, 


! / / 1 A X , 4 ( 3HNEW , 7 X ) , 7 X , 4 ( 3H0LD , 7 X ) / , 


!7H NODE , 2(4X, 42HU-VEL0CrrY PRESSURE V-VELOCITY TENPARATURE) ) 
DO 20 IP0IN=1, MPOIN 
IODFM=NODFM(IPOIN) 


IADFM"NADFM( IPOIN) 

IF ( lODFM.EQ. 'NDOFM) GOTO 10 

WR I T E ( 35 , 20 1 0 ) I PO I N , ( VARB 1 ( I ADEN ^ JODFN ~ 1 ) , JODFM= 1 , I ODFM ) , 
! (VARB2(IADFM+J0DFM-1), JODFM^l, I ODFM) 

DO J0DFM=1, lODFM 
JK-^JODFM 

IF (JK.GT. 1)JK"JK+1 

KK = IADFMi-JODFN-l 

WRITE(40, «) IPOIN, JK, VARBl (KK) 

IF(JK.E0.4)THEN 

VARBX ( KK ) = VARB 1 ( KK ) «T REF 

WRITE(88, ») IPOIN, VARBX(KK) 
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E^4D I F 
FMDDO 
60 10 20 
10 CONTINUE 

WRITE (35, 2020} IPO IN, (VARBl ( I ADFM+JODFN-1 ) , J0DFM=1 , lODFMI , 

! (VARB2(IADFM+vJ0DFM-l}, J0DFM=1, I0DFM1 
DO J0DFM=^1, lODFN 
JK-^-^JODFM 

I I ADFM+JODFN-1 

WR I TE ( 40 , w ) I PO I N , JK , VARB 1 ( KK ) 

IF(JK.EQ. 4) THEN 
VARBX(KK.)=VARB1 (KK)sTREF 
WR I TE ( 88 , w ) I PO I N , VARBX ( KK ) 

ENDIF 

ENDDO 

20 CONTINUE 

2010 FORMAT (IA,4X,2(E10.3, 10X,2E10.3,6X) ) 

2020 F0RMAT(I6,4X,?(4E10.3,6X)) 

RETURN 

END 

C 

SUBROUTINE TOLREL( I ITER, NTOTV, NCONV, NTOTV, RELAX , TOLER, 

! VARBl , VARB2, NADFM, NODFM, MPOIN) 

IMPLICIT REAL«8(A-H,0-Z) 

D I MENS I ON VARBl ( 2000 ) , VARB2 ( 2000 ) , NADFM ( 386 ) , NODFM ( 386 ) 

C TO CHECK IF SOLN. HAS CONVERGED 
ND0FM--4 

DO 200 IPniN=l , MPOIN 
lODFM^NODFMdPOIN) 

IADFM=NADFM(IPOIN) 

IF(IODFM.EQ.NDOFM) GOTO 108 
GOTO 208 
108 CONTINUE 
208 CONTINUE 
CANLA=0. 0 
NCONV-'l 
UMAX=VARB1 (1) 

PMAX^VARBl (2) 

VMAX-VARBl (3) 

TMAX=VARB1 (4) 

DO 10 IP0IN==1, MPOIN 
lODFM-MODFM(IPOIN) 

IADFM=NADFM( IPOIN) 

IFdODFM.EO.S) GOTO 16 

IF (DABS (VARBl ( lADFM) ) . GT .DABS (UMAX ) .1 UMAX=VARB1 (lADFM) 
IF(DABS (VARBl ( IADFM+1 ) ) . GT. DABS(PMAX ) ) PMAX=VARB1 (lADFM+l ) 
IF (DABS (VARBl ( IADFM+2) ) . GT. I)ABS(VMAX ) ) VMAX-VARBl ( IADFM+2.') 
IF (DABS (VARBl ( IADFM+3.) ) . GT. DABS(TMAX) ) TMAX=VARB1 dADFM+3) 
GOTO 100 

16 IF (DABS (VARBl (lADFM) ). GT. DABS (UMAX ) ) UMAX=VARB1 (lADFM) 

IF( DABS (VAREil ( I ADFM+l ) > . GT. DAB3(VMAX ) ) VMAX=VARB1 ( IADFM+1 ) 
IF (DAPS ( VARFil dADFM+2) ) . GT. DABS(TMAX) ) TMAX--VARB1 dADFM+2) 
100 CONTINUE 

'10 CONTINUE ' 
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UCUT^O. IsUMAX 
PCUT=0. IwPMAX 
VCUT=0. IwVMAX 
TCUT=-0. 1«TMAX 
DO 18 IPOIN=l,MPOIN 
IODFM=NODFM(IPOIN) 

IADFM=NADFMC IPOIN) 

ITOTU=IADFM 

IF( lODFM.EQ. 3) GOTO 19 
ITOTP^-^IADFMh 1 
ITOTV=ITOTU+2 
itott=itotu+3 
GOTO 34 

19 ITOTV^-TTCrfUfl 
ITOTT-.:-^rTOTU<-2 

34 I F (. DABS (. VARP 1 ( ITOTU ) ) , LT . DABS (UCUT ) ) GOTO 60 
IF(DABS(VARB1 (IT0TU>3.LT.0.01) GOTO 60 
CANGF"DABR( (VARBl (ITOTU) -VARB2( ITOTU) ) /VARBl (ITOTU) ) 
I F ( CANGE , GT . TOLER ) NCaNV=0 
IF (CANLA.GT. CANGE) GOTO 60 
I CANI.A-^CANGE , 

L..TOLA= ITOTU 
60 CONTINUE 

IF( JODFM. EO. 3) GOTO 70 

IF(DABS(VARB1 (ITOTP)).LT.DABS(PCUT)) GOTO 70 
I F ( DABS ( V ARB 1 ( I T OTP ) ) . LT . 0 . 0 1 ) GOTO 70 
CANGE^DABS( (VARBl ( ITOTP) -VARB2 ( ITOTP )) /VARBl ( ITOTP) ) 
I F ( CANGE . GT . TOLER ) NCONV-0 
IF(CANLA.GT. CANGE) GOTO 70 
CANLA=CANGE 
LTOLA= ITOTP 
70 CONTINUE 

I F ( DAPS ( V ARB 1 ( I TOT V ) ) . LT . DABS ( VCUT ) ) GOTO 80 
IF ( DABS (VARBl (I TOTV) ) .LT. 0,01) GOTO 80 
CANGE=DABS( (VARBl ( I TOTV) -VARB2( I TOTV) ) /VARBl ( I TOTV) ) 
I F ( CANGE . GT . TOLER ) NC0NV=0 
, IF (CANLA.GT. CANGE) GOTO 80 
CANLA=CANGE 
LTOLA=ITOTV 
80 CONTINUE 

IF(DABS(VARB1 ( ITOTT) ) . LT.DAB5(TCUT) ) GOTO IS 
IF(DABS(VARB1 (ITOTT) ) . LT. 0. 01 ) GOTO 18 

CANGE="“DABS( (VARBl ( ITOTT) -VARB2 ( ITOTT ) ) /VARBl (ITOTT) ) 
I F ( CANGE . GT . TOLER ) NC0NV=0 
I F ( CANLA . GT . CANGE ) GOTO 1 8 
CANLA-^CANGE 
LTOLA= ITOTT 
18 CONTINUE 

WR I TE ( , 2000 ) LTDLA , CANL A 
WRITE(35, 2000)LT0LA , CANLA 

2000 FORMAT ( /37X, ■' LARGEST CHANGE OCCUR AT DOF NO ',15,/, 
!10X,' CHANGE =',F10.4) 

IF(NCONV.Ea.O) GOTO 20 
WRrrE(M,2010) IITER 
WR I TE ( 35 ,2010)1 ITER 


2010 FORMAT (//^5H SOLUTION HAS CONVERGED TO REQUIRED TOLERANCE/, 
!3X, 'INSSX, 15, 3X, 'ITERATIONS') 

RETURN 
20 CONTINUE 

IF (I ITER. EQ. 1.) GOTO 40 
C 

C IJIELAX VAR. FOR NEXT ITERATIONS 
DO 222 IPOIFT^-l ,MPOIN 
IODFM=^^^NODFM(IPOIN) 

IADFM=NADFM( IPOIN) 

DO 21 JODFN=a, lODFM 
REL=RELAX 

I F ( C I ODFM . EO . 4 ) . AND . ( JODFM . EQ . 2 ) ) REL= 1 . O 
TEMP Y 1 =VARB 1 ( I ADFM+ JODFM - 1 .) 

TEMPY2=VARB2<: lADFM+JODFM- 1 ) 

VARB2(IADFM+J0DFM-l.)=TEMPY2+RELw(.TEMPYl -TEMPY2) 

,21 CONTIN\)E 
222 CONTINUE 


40 CONTINUE 

DO 50 IT0TV=1,NT0TV 
VARB2( I rOTV.)-VARBl ( ITOTV.) 

50 CONTINUE 
225 CONTINUE 
RETURN 
END 

C 

SUBROUTINE FRONTS (.lAXSY, I ITER, MELEM, MFRON, MPOIN, MTOTV, NBCON 
! LEM , NEVAB , NNODL , NNODP , NPO IN , NTOTV , XFORC, YFQRC , NGAUS , NESPBC , I 
!M, BOUDV, LBOUD, LHEDV, PNDRM, EQFIHS) 

IMPLICIT REAL»8CA-H,0-Z) 

C EXTTERNAL SUBROUTINES 

C MATRIX 

DIMENSION LBOUD (2000) , LHEDV ( 119), PNORMd 19) , BOUDV (2000) 

D I MENS I ON EORHS ( 2000 ) 

C0MM0N/ARC2/I.. NnDS( 1 1 1 , B) , CnORD(3B6, 2) , NADFM(386) , N0DFM(3B6) 
C0MM0N/ARC3/VARB1 (2000) , VARB2(2000) 

D I MENS I ON LOCEL ( 28 ) , NDEST ( 28 ) , FLUMX ( 28 , 28 ) , GFLUM (119,119) 
OPEN (UNI T~25, FORM= •'UNFORMATTED ' , FILE= 'ASDT3. OUT-' ) 
IF(KELEM.E.Q.NELEM) GOTO 67 
67 CONTINUE 

IFdITER.GT. DGOTO 40 

C ON FIRST ITERATION ONLY FIND LAST APPEARNCE OF EACH NODE 
DO 30 IFOIN^-1 ,NPGIN 
I ASTF""0 

DO 20 IELFM^-1,MFLEM 
DO 10 IN0DP=1,NN0DP 

IF(LMODS(IELEM, INODP) .NE. IPOIN) GO TO 10 

laste^v,= i[.:lem 

1 LASTN-INODP , 

GOTO 20 
10 CONTINUE 


,NE • 
GFLU 
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?o cnr-iriNUE 

I M( tDB ( I. ABIE , LAST N ) ~ - 1 I N 
30 CONTINUE 
40 CONTINUE 

REWIND 23 
NCniT^'^MFRON- NEVAB 
NFR0N--0 

DO 50 ir FU.:)M=-3. ,MF RGN 
DO 50 JFR0N=1,MFR0N 
or I . UM (. I FRON , JFRON ) =0 . 0 
50 CONTINUE 
KELEM--^0 

C START ASSEMBLY BY FORMING ELEMENTAL MATRICES 
, 60 CONTINUE 

KFI. FM--KELEMM ' 

CALI. MAT FU X ( I AXSY, MELEM, EORNS, FLUMX , MPOIN 
! , MTOTV, NEVAB, NNODL, NNODP , XFORC, YFORC 
! , KELEM, NELEM, NGAUS, I ITER ) 

C IF(KELEM.EQ.NELEM)GOTO 3770 

C 3770 RETURN 
KEVAB=0 

C CREATE GLOBAL DOF ARRAY FOR EACH LOCAL ELEMENT DOF 
DO 70 INODP^lrNNODP 
KPOI N=LNODS ( KELEM, I NODE .) 

I ADFM=NADFM( lABS(KPOIN) .) 

LODF- M^='^NODFM ( I ABS (KPOI N ) ) 

DO 70 lODFM^l ,LODFM 

KEVAB=KEVAB+1 

LOCEL ( K E V AB .1 = I ADFM+ 1 0 DFM - 1 

I F ( KP 0 1 N . LT . 0 ) LOCEL ( K EVAB ) = -LOCEL ( K E V AB .) 

' 70 CONTINUE 

C FIT EACH DOF INTO THE FRONT WIDTH EXTENDING IF NECESARY 
DO 120 lEVAB-^i ,NEVAB 
KTOTV-=LDCEL(IEVAB) 

TF(NFR0N,Ea.0,>G0T0 90 
DO 80 IFRON^lrNFRON 
KFRON-IFRON 

IF( IABB(KTOTV.I .EQ. I ABS (LHEDV (KFRON.) ) .I GOTO 110 
80 CONTINUE 
90 CONTINUE 

NFR0N--’NFR0N+1 

IF (NFRON.LE.MFRON)GOTO 100 
WRITE(»,2000,) 

2000 FORMAT!. //3X, ■•'PROGRAM HALTED FRONT WIDTH IS TOO SMALL'') 
STOP 

100 CONTINUE 

NDEST ( I EVAB) =NFRON 
LHEDV (NFRON) -KTOTV 
GOTO 120 
110 CONTINUE 

NDEST ( I E VAB ) =KFRON 
LHEDV (KFRON)=KTOTV 
120 CONTINUE 

C ASSEMBLE NEW ELEMENNT IN TO GRAND FLUID MATRIX 
DO 1 30 I EVAB= 1 , NE VAB 
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f -J N 


I 


I.FROM’==NDEST( IEVA8) 

I)n 130 .TFVATV-n ,Nr VAB 
JFROM-^JDESl'c.JFVAB) 

GFLUM ( JFRON, IFRON)=BFLUM(aFRON, ir RON)-! FLUMX( JEVAB, lEVAB) 
130 CONTINUE 

I F ( MF ROM . LT . NCR I T . AND . KELEN . LT . NFIEM ) GOTO 60 
140 CONTINUE 
NF SUM- 0 


PIVOT-0.0 

C CHECK LAST APPEARANCE OF EACH DOF PROCESS BOON. CODNS. 

DO 1 /O I FROM- 1 , NFRON 
ir(I..HEDV(IFRON) .GE.OIGOTO 170 
NFSUM-=1 

I F ( LBOUD ( I ABS ( LHEDV ( I FROM ) ) ) . NE . 1 ) GOT 0 160 
KTOTV^=IABS<,LHEDV( IFRON) ) 

LBOUD (KT0TU)“-1 

EORHS ( K rOTV T ^=BOUDV ( K TOTV ) 

DO 130 L F ROM™- 1 , NFRON 
GFLUM ( I FROM , LFRON ) =0 . 0 
150 CONTINUE 

GF I ,UM ( I FROM , I FROM ) = 1 . 0 
160 CONTINUE 

C SEARCH FOR LARGEST PIVOTAL VAL 
P IVOG“-GFLUM( IFRON, I FRONT 
IFIDABSIPIVOG) .LT.DABS(PIVOT) TGOTO 170 
PIVOr'=PIVOG 
LPIVT"IFRON 
170 CONTINUE 

I F (. NF SUM . EQ . 0 ) GOTO 60 
K TOTV" I ABS ( LHEDV ( 1_P I VT ) ) 

I F ( DABS (. P I VOT ) . GT , 1 . OE -25 ) GOTO 1 80 
C I F ( DABS ( P I VOT . GE . . OE+00 ) GOTO 1 80 

WR I TE ( w , 20 1 0 ) KTOTV , P I VOT 

2010 FORMAT! //'PROGRAM HALTED ILL -CONDITIONING, '// 'D. 0. FREEDOM’' 
! , 14, / -'PIVOT VALUE-', E9. 2) 

STOP 

180 CONTINUE 

Cw >>: •o-sf-ft! -'j; -o! NORMALISE. F'lVOTAL E.QN 
DO 190 I FROM" 1 , NFRON 

P NORM (. I FROM ) "OF I...UM ' L.P I VT , I T RON .> / P I VO T 


C 


190 CONTINUE 

RHS I D^’T:0F^,HS (.K TO tv .T /P I VOT 


EORHS ( K TO TV ) ^-^RI-IS I D 


«^<wELIMINA"f ION OF PIVOTAL EON. REDUCING FRONT WIDTH 
, IF(LPIVT.EQ. l.TGDTO 250 , 

DO 24 0 IF R^M=^•^1, ,LPIVT-1 


FACOR"GFLUM( IFRON,LPIVT) 

I F ( DABS ( FACOR ) . LT . 1 . OE -30 ) GOTO 2 1 0 

DO 200 JFR0N=1,LPIVT-1 

GFL UM ( I FROM , JFRON ,1 =6FLUM ( I FRON , JFRON .1 


-FACDR^^PNORMf.JFRON 


) 


00 

10 


CONTINUE 
COMT I NUE 

I F ( LP I VT . E Q . NFRON .> GOTO 230 
DO 220 JFRON ’M_P I VT+1, NFRON 
GFLUM (. I FRON , JFRON - 1 ) "GFLUM <. I FRON , JFRON ) 


-FACORwPNORM ! JFRON) 
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220 CONTINUE 
230 CONI T NUE 

n OT I ABS ( L HEDV Cl FRON > ) 

EQRHS C I TOTV ) ^EQRHS C I TOTV ) -FACOR^^RHS T D 
240 CONTINUE 
250 CONTINUE 

IFCLPIVT.EO.NFRON)BOTn 300 
■ DO 270 IFR0N=LPIVT+1 ,NFRON 
FACOR=^ RFL.UM C I FRON , LP I VT ) 

IF CLP I VT. HO. 1 )GOTa 270 
DO 260 JFR0N=1,LPIVT-1 

or I UMC IFRON-l , JFRaN)=GFLUMC IFRDN, JFRDN) -FACORsPNDRMC JFRON). 

260 CONTINUE 
270 CONTINUE 

DO 2R0 JFR0N=LPIVT+1,NFR0N 

GF L, UN C I FRON - 1 , JFRON - 1 ) =GFLUM ( I FRON , JFRDN ) -F ACORkPNDRM ( JFRON ) 
280 CONTINUE 

I T07V“ I ABS ( LHEDV C I FRON ) ) 

EORHSC ITOrV)=EORHS(ITOTV) -FACORwRHSID 
290 CONTINUE 
300 CONTINUE 

C?^k«WRITE NONFIXED pivotal EON. 

IFCLBOUDCKTOTV) .NE.OIGOTO 310 

WR I TE C 25 ) NFRON , LP I VT , C LHEDV C I FRON ) , PNORM C I F RON ) , I FRON= 1 , NFRDN ) 
310 CONTINUE 

DO 320 IFRDN=^1,NFR0N 
GFLUMC I FRON, NFRON )=0. 0 
GFLUMC NFRON, I FRON) =0.0 
320 CONTINUE 

IF CLPIVT.EGl. NFRON) GOTO 340 
DO 330 IFR0N=LPIVT,NFR0N-1 
LHEDV C I F RON ) =L.HEDV ( I FRON+ 1 ) 

330 CONTINUE 
340 CONTINUE 

NFR0N=NFR0N-1 

C5<«-»-« ASSEMBLE ELIMINATE OR BACK SUBSTITUTE 
I F C NFRDN . GT . NCR I T) GOTO 1 40 
I F C K ELEM . LT . NELEM ) GOTO 60 
I IF (NFRON. GT. 0) GOTO 140 , 

SUBSTITUTION 
DO 350 IT0TV=1,NT0TV 
VARB KIT OTV ) =BOUDV ( I TOTV ) 

LBOUDC IT0TV) = --LB0UD(IT0TV) 

■350 CONTINUE 

DO 370 IT0TV=1,NT0TV-NBC0N ' 

BACKSPACE 25 

READ (2.5)NFR0N, LP I VT, (LHEDV ( IFRDN) , PNORM ( I FRON) , IFR0N=1 , NFRON) 
KTOTV=IABS(LHEnV(LPIVT) ) 

TEMPR=0.0 
PNORM(LPIVT)-=0.0 
DO 360 IFR0N=1,NFR0N 

TEMPR=’ T EMPR -PNORM ( I FRON) sVARBl ( I ABS (LHEDV ( IFRDN) ) ) 

360 CONTINUE 

VARB 1 (KTOTV ) =EDRHS (KTOTV) +TEMPR 
, BACKSPACE 25 
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370 CONTINUE 
RETURN 
END 




C 

SUBROUT I NE NATR I X Cl AXSY, MELEM, EQRHB , FLUMX , 

! MPOIN, MTOTV, NEVAB, NNODL, NNODP , XFORC, 

! YFORC, HELEN, NE LEM, NBAUS, I ITER) 

IMPLICIT REALx8CA-H,a-Z) 

C«*^«INSERTED THE DUMMY VARIABLE "NSIDE" APOVEhwwwwwwwwsw 

DIMENSION CARTLC2,4) ,CARTP(2,8) ,ERHSUC8) ,ERHSV(8) ,FLUMXC2a,2a) 
DIMENSION ERHSTC8),VISCC9),SHAPXC8) 

D I MENS I ON SHAPL C 4 ) , SHAPP C 8 ) , AREAW C 9 ) , DERX (2,8), DER I V ( 2 , 8 ) 
C0MM0N/ARC2/LN0DS( 1 11,8), CQ0RD(38A, 2) , MADE MC 386) , N0DFM(386) 
DIMENSION EaRHS(2000) 

C0MMDN/AREA3/CARPB(2,72) , SNAP G( 72) , CARLS (2, 36) , SHALB(36) 
C0MM0N/ARC1/P0SGP(3) ,WEIGP(3) 

C0MM0N/ARC3/VARB1 (2000) , VARB2(2000) 

C C0MMaN/ARC4/AReAW(9) 

COMMON/ ARC 1 0/CUTVEL , DCUT , RAKE 

COMMON / ARC 1 2 / FR I FAC , P I , EN , COS A , S I NA , S I N2 A 

CaMM0N/ARC13/RH'0, MBOUEL 

COMMON /ARC 14/1 BEL ( 60 ) , I SURF ( 60 ) , I S I DE ( 60 ) 

COMMON/ ARC35/EPSN ( 1 1 1 , 9) , XB ( 1 1 1 , 9) , YG ( 1 1 1 , 9) , VISCP (111,9) 

! , DARE ASCII 1,9) 

COMMON/ ARC99 1 /TREE , PECLET , CONRAT , 1 FST AR, H57 AR , 
!BSTAR,SIGYO,SPHTO,THKO, SIGMA K.K.,-»Ar, 

CALL DRIVES (COORD, LNODS, MELEM, MPOIN, NBAUS, NELEM, NNODL, NNODF , 

! HELEN, AREAW) 

C 

C«»ww INITIALISE ARRAYS 

' DO 10 IN0DP=1,NN0DP / 

ERHSU(INODP):='^0.0 

ERHSV(INODP)=0.0 

ERHST(INODP)=0.0 


10 CONTINUE 

DO 20 IEVAB=1,NEVAB 
DO 20 JEVAB=1,NEVAB 
FLUMX ( lEVAB, JEVAB)=0. 0 
20 CONTINUE 


C 

C LOOP TO CARRY OUT GUASS INTEGRA f IONS 


C 


IGAUS=0 

DO 999 IB=1,3' 
DO 999 JG=1,3 
DO 100 IGAUS=1,9 
IGAUS=IGAUS+1 
XEQIV-“POSGPC.IG) 


YE0IV=P0SGP(J6) 

CALL SHAPEBCDERIV, SHAPX , XEOI V, YEOIV) 


TEMPR.=0. 


XG(KEi:LEM, IGAUS)=0. 
YGCKELEM, IGAUS)=0. 


I 
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DO 222 INODP^=l,a 
K P 0 1 N-- 1 ABS ( I NODS ( K EiiL ,EN , I NGDP ) ) 

ITOTU ^NADFN(KPOrN) ' ' 

I I f ri' T I 'i cn u t nc:)Df: m ( k ? n i n ) 1 
Tt£MPR^^^-'.TEMPR'4 VARB2( HOT D wSHAPX f INQDP) 

C FINDIND GLOBAL POS. OF G.PTS 

XG(Kt;l Eli, IGAUS)=XG(.KELEM, IGAUS) + (SHAPX ( INODP ) sCDDRD f KPOIIM, 1) > 
! sDCUT:^1000. 

YG ( K E.L, EM , I GAOS ) Y G ( K ELEM , I GAOS ) + ( SHAP X ( I NDDP ) sCOORD ( K PO I N , 2 ) ) 
IwDCUTwlOOO. 

222 CONTIfJUE 
C 

I CALL PROP (TEMPR,SI6STAR,CPSTAR,THKSTAR,SIBMAY) 

C 

DO KINP=0,fJ 

KPfllN-IABSd NGnS(KEl.EM,KINP) ) 

ITOTU~NADFM(KPOIN) 

I T OT T= I TO’I Of NODF M ( KPO I N ) - 1 

KIGS^^HGAUS 1 

I)AREA=^AREAWa<IGS) 

DAREAG ( KELEM , K I GS ) ^^-DAREA 
DO 30 :iNODP=..l,IMNni}P 

SHAPl^ ( I NDDP ) -SI lAPG ( NNODP ( I GAUS - 1 ) + 1 NODP ) 

DO 30 ID1ME=H2 

OART P ( I DIME, INODP ) --^CARPG ( IDIME, NNODP« ( IGAUS-1 ) + INDDP ) 

30 CONTINUE 

DO 40 INODL=sl,NNODL 

SHAPE ( I NODE ) -SHAEG ( NNODl . w ( I GAUS - 1 ) H NODE ) 

DO 40 I DIME- 1,2 

C;',ARTL n D I ME , I NODI. ) -C ARLG (I D I ME , NNODEw ( I GAUS - 1) + 1 NODE ) 

40 CONTINUE 

IEEEM--KELEM 
UVEEY-O.O 
RADUS-0. 0 
WEEY-O.O 

Can EVAE. RAD.S.:PREV. UEES AT GUASB POINTS 
DUD X -0.0 
DUDY-0.0 
DVDX=0.0 
DVDY-0.0 
C 

DO 30 I NODP ==1, NNODP 

KPO I H~ 1 ABS ( ENODS ( KEEEM , I NODP ) ) 

rrOTU~NADFM(KPOIN) 

I TOTV-I TOTUH-NODFM (KPQIN) -2 
SHAPE-SHAPPdNODP) 

UVEEY^=UVEEY + VARB2( ITOTU.) wSHAPE 
HADUS=T^ADUSh COORD (KPOIN, 2) itSHAPE 
VV£EY^'=VVEEY+VARB2 ( ITOTV.I sSHAPE 
CARX I^CARTP ( 1 , INODP .) 

CARYr^CARTP(2, INODP 
UVEL-VAflbDnTOTU) 

VVEE'--VARB2( rrnru) 

DUDX^DUDX+CARX HUVEL, 

DUDY -DUD Y * CARY I s^UVEL 
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DVDX ^=DVDX iCARX I ^ WEL 
DVDY^^ DVDY-i- CARY I iiVVEl.. 


C 

50 CONTINUE 


C 

EPB^T)BQRTt2. 0s(DUDXii:Ji2+0.5iit (DUDY+DVDX)s«2.HDVDYicH2) ) 
C 

I F ( EP S , I .T -0.0001) EP B=-0 .0001 
EP BN ( K EL EN , K I G3 ) =^EP 3 


C 

V I SC ( K I GS ) = (. S I GST AR+BST ARj^EPSs^ (1 . 0/EN ) ) / ( 1 . 732ftEPS ) 

V I SCP ( K ELEN , K I GS ) ^ V I SC ( K I GS ) 

C 

IFdAXSY.EQ. 1 )DAREA=DAREAteRADUS 
Ciiws PUT HEAT GEN. TERM/OR BODY FORCES IN TO LOCAL RHS VECTOR 
C 

DO 60 INni)P~U,NNaDP 

ERNST t INODP ) =^T: RHST ( INQDP ) +SHAPP ( INODP ) «EPS«PECLET»DAREA 
! wEPS:<iVISC(,KIG3.) 

60 CONTINUE 
' DO 90 IC0N1=1,4 

DO 90 IC0N2=-1,2 
I ROWU~ ; I COM 1 - 1 ) ■<t 7+-4s I C0N2 -3 
I ROWV-- 1 RUWUt-3 •- 1 CGN2 
I ROWT=^ I R0WU <-4 - 1 C0N2 
I NODP" 2« ( I qON 1 1 ) + 1 C0N2 
SNAP I "SHARP ( INODP) 

CAHXI=-CARTPC1, INODP) 

CARY I --C ARTP (2,1 NODP ) 

IF (IC(;)N2. EG. 1 ) IR0UP=IRDWU+1 
DO 80 JCON 1-1,4 
JCOLP = ( JCON 1 1 ) k 7 -t'2 
K GAL I = ( K I GS - 1 ) wNNODL f JCON 1 


C 

C«MW PUT PR. TERM IN MOM. EQN. ,x, r- , ,,rvA. T^ 

FLUMX ( IROWU, JCOLP ) --FLUMX ( IROWU, JCOLP ) -CARX IwBHALG (KGALI ) 

* '"''^FLUMX (IROWV, JCOLP) -FLUMX (IROWV, JCOLP) -CARYIj^BHALeCKGALI ) 
I »DAREA 

DO 80 JC0N2-1,2 
JCOLU~( JCONl -1 ) W7+4-4JC0N2 -3 
JC0LV-JC0LU-+3”JC0N2 
JCOLT" JCOLU+4 - JC0N2 
JN0DP-2«( JCONl -1 )+JC0N2 
SNAP J=SHAPP ( JNODP ) 

CARX J=CARTP ( 1 , JNODP ) 

CARYJ--CARTP ( 2 , JNODP ) 

P^'T WFFU ^ CONVE TERMS IN 

rnmc-iaml j+wei.y^cary j ) (shap i .dareai > /si eha 

cfumT HOUU irOLU) "FLUMX ( IROWU, JCOLU) +DIFFU+CONVC 
m ™ CARX I KMRXa+2 . o"cARY I «CARY J ) I SC ( K 1 8S ) *DAREA 
c, imvT mnwJ SrolO ) UMX 1 1 RDWV , JCOLV ) 4-01 FFU+CONVC 
fumaSu:™LV;"cUl*CABXT*VISCIKIGS,.I.ARE^^ 

! +FLUMX ( I ROWU , J COLV ) 
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F L UMX ( I ROWV , JCOLU ) =FLUIiX (. I nOWV , JCOLU ) +CARX I i^CARY I SC ( K I GS ) 

! i^DARFA 

cmww form energy eon 

CwswwPUT COND. S.:CONV. TERMS IN ENERGY EON 

EMCOND= (. CARX I sCARX J+CARY I «CARY J ) sUARE As T HKST AR 
ENCONV= (UVELYsCARXJ+VVELYsCARYJ) sSHAP I sDAREAsPECLET 
IsCPSTAR 

FLUMX (IROWT, JCOLT)=FLUMX ( IRDWT, JCOLT)+ENCOND+ENCONY 
IF(IC0N2.EQ.2) GOTO 70 
psss FORM CONTINUITY EON 

FO....UMX (. IROWP , JCOLU) =FLUMX ( IROWP, JCOLU) +SHAPL( ICONl ) sDAREAsCARXJ 
FLUMX ( IROWP , JCOLU) =FLUMX ( I ROWP , JCOLV) +SHAPL( ICONl ) sDAREAsCARYJ 
FLUMX ( I ROWP , JCOLP ) =0 . 0 
FLUMX ( I ROWP , JCOLT ) =0 . 0 
70 CONTINUE 
80 CONTINUE 
90 CONTINUE 
ENDDO 

999 CONTINUE 

CswssssssssssADD LOCAL RHS VECTOR TO GLOBAL AFIRAY 
DO 110 INaDP=l,NNODP 
KINP=INODP 

KPOIN=aABStLNQDSnELEM,KINP) ) 
rraTU=NADFM(.KPOIN) 

I TOT V= I TOTU+NODFM f. KPO I N ) -2 
ITOTT^ITOTU+NODFMCKPOIN) -1 
EORNS ( I TOTU ) “EORHS ( I TOTU) 4ERHSU ( INODP ) 

EQRHSI ITOTV)=EQRHStrrOTV) t-ERHSV(INQDP) 

EORHS ( I TOTT ) =--EQRHS ( I TOTT ) +ERHST ( I NODP ) 
no CONTINUE 

DO 120 I=1,MB0UEL 
JK=I 

IFIKELEM.NE. IBELIJK) ) GOTO 120 
KBEL=IBEL(JK) 

KSURF=ISURF (JK) 

C 

I F (. ( KSURF . EQ .-5 ) . OR . ( KSURF . EQ . 6 ) ) THEN 

DO 1001 IC0N1=1,4 

DO 1001 IC0N2=n2 

I ROWV= I ROWU+3 - I C0N2 

I N0DP=2s ( I CON 1 - 1 ) + 1 C0N2 

IF( IN0DP.LG.3) THEN 

DO 1002 JC0N1=1,4 

JCOLP- ( JCON 1 -l.)»7+2 

DO 1003 JC0N2=1,2 

JC0LU=(JC0N1 - l)s7+4sJC0N2-3 

JCOLV^^ JCOL.U+3 “ JC0N2 

JCOLT" JCOLU+4 -JC0N2 

JN0DP-2S ( JCON 1 " 1 ) + JC0N2 

fUJJMX ( IROWV, JCOLU) “SINAsFLUMX ( IROWU, JCOLU) +- 
! COSAsFL UMX (. I ROWV , JCOLU ) 

FLUMX ( IROWV, JCOLV)=SINAsFLUMX ( IROWU, JCOLV)+ 

! COBAsFLUMX ( IROWV, JCOLV) 

FLUMX t IROWV, JCOLT) ^SINAsFLUMX ( IROWU, JCOLD f 
! COSAsFLUMX (, I ROWV, JCOLT ) 

I 
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u u o u o u u 


/ 


1003 COMTIMUE 

F I.UMX ( IROWV, JCOLP >“BIMA«rLUMX ( TROWU, Jf'Cll P) + 
ICOSAy^FLUMX ( IROWV, JCOLP) 

1002 CONTINUE 
ENDIF 

1001 CONTINUE 
END IF 

IP ( (KSUF1F.LE.3).DR. IKBURF . FO. 7 ) lOOTD 170 
KSIDE^=ISIDE(JK) 

CAI .1, BOUCOM ( K BEL., K surf , KSIDE, FLUMX , EORFIS , VARB2, MELEM, 
! MPO I N , NNODP , MTOTV , 1 1 TER .) 

120 CONTINUE 
RETURN 
END 


C W ■)’: W M W W W W W W W W W W W W H M i’: -K •t': W W W 'W W ii< W W !n- W W W W W W i!*: W W W 


I 


SUBROUT I NE DR I VES ( COORD , LNODS , MELEM , MPO I N , NELEM , NGAUS , 
! NNODl. , NNODP , I ELEN, ARE AW, 1 
I NPL I. C I T REALMS ( A -H , 0 -Z .F 


W W W W W K' 1^' W W W W W M W iof S W -ot W S W W 


EXT E'UBROUTINES 


DJACOB 

D I MENS I ON COORD ( 386 , 2 ) , L NODS ( 1 1 1 , El .1 
D I NENS I ON DER I V ( 2 , 8 ) , DJAC I ( 2 , 2 .) , D JACK ( 2 , 2 , SHAPE ( 8 ) 

D I NENS I ON C ART P (2,8), ARE AW ( 9 ) 

COMNON/ARC 1 /POSGP (3 ) , WE I GP ( 3) 

CONNON/ AREA3/CARPG (2,72) , BHAPG (72) , CARL.G(2, 36) , SHALG(36) 
C0NN0N/AREA4/AREAW(9) 

BET UP POSITIONS ?, WIS FOR 3 POINT GUABS RULE 
POSGP ( 1 ) ^-0. 7745966692 
POSGP (2) =-0.0 
P0SGP(3)=-P03GP(1) 

WE I GP ( 1 ) ^0 . 5555555556 
WE I QP ( 2 ) =0 . 8888888889 
WEIGP(3)=WEIGP(1) 

C CALCULATE SHAPE FUNCTIONS S-: DERIVATIVES FOR ELENENTS 
LGAUS=^0 
NGAUS=3 

DO 50 IGAUS=1,NGAUS 
DO 50 JGAUS-1 , NGAUS 
LGAUS=L.GAUS+1 
XEQIV=POSGP(IGAUS) 


YEQIV==f'GSGP(JGAUB) 

C USE GAUS POSITIONS TO CALCULATE LOCAL VALUES 
CAl.L SHAPES (DERIV, SHAPE, XEQIV, YEQIV) 

C SET UP JACOBIAN MATRIX & INVERSE 

CAIJ... DJACOB(DERIV,DETJB,DJACI,DJACK, lELEM, 
! MEI ..EM , MPO I N , NNODP , LGAUS ) 

C CALCULATE GLOBAL.. DER t. AREAwGAUS WTS 
DO 10 IDIME=lr2 
DO 10 IN0DP==1,NN0DP 
CAR TP ( I D I ME , I NODP ) =0 ; 0 


DO 10 JD I ME" 1,2 

CART P ( I D I ME , I NODP ) =T; ART P ( I D I ME , 


I NODP ) +D J AC I ( I D I ME , JD I ME ) 
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- !?^DFRIVf.JDIME, INODP) 

10 CONTINUE 

AREAW (. LBAUS ) DET JBwWE I GP ( I GAUB ) wWE I GP < JGAUB ^ 
C PUT SHAPE FNS ?•: pERIVB IN ELEMENT MX 
DO 20 IN0DP=1,NN0DP 
KARA - (LGAUS-1 ) wNNODPHNODP 
BH AP G (. K AP A ) AP E ( I NODP ) 

! DO 20 IDIME=1,2 

CARPG ( I D I ME , K AP A ) =CARTP ( I DI ME , I NODP ) 

20 CONTINUE 

C USE GAUSS POSITIONS TO CALCULATE LOCAL VALUES 
CALL SHAPE4(.DER IV, SHAPE, XEQIV, YEQIV) 

C CALCULATE GLOBAL DERIVATIVES FOR LINEAR FUNCTIONS 
DfJ 30 rDIMF = l,2 
DO 30 IMaDL=l,MNaDL 
CAR T P ( I D I ME , I NODE) =0 . 0 
DO 30 JDIME=a,2 

CART P f. 1 1.) I ME, , I NQDL ) “CARTP ( I DIME , I NODL ) + 

! D JAC I U D I ME , JD I ME ) wDER I V (. JD I ME , I NODL I 
30 Cf,)MTINUE 

C PUT SHAPE: FNS P. DERI VS IN ELEMENT MX 
DO 40 IN0DL:=1,NN0DL 
KBAL. I ™ ( LGAUS -■ 1 ) ttNN0DL+ 1 NODL 
SHALGTKGALI ) -SHAPE ( I NODL I 
40 CONTINUF.:; 

50 CONTINUE 
RETURN 
END 

Q ^ ^ ^ ^ ^ ^ ittr ior W ior iBr itj: -ct ice « -oj to? W ^0: ior K ioJ W ior W 

c 

BUBROUT I ME SHAPES ( DER I V , SHAPE , XEQ I V , YEO I V ) 
IMPLICIT REALMS ( A -H,0-Z) 

DIMENSION DERIV(.2,8,') ,BHAPE(8) 

C PARABOLIC ELEMENT ANTICLOCKWISE NODE NUMBERING 
X=XEQIV 
Y==YEQIV 
XY-X»Y 
XX-^X;«X 
YY“Y«Y 
XXY^XXwY 
XYY=-^^X?‘YY 
X2""Xw2. 

Y2-Y«2. 

XY2-XYW2. 

shape; ( 1 ) = I “1 . ■+ XY+XX+YY“XXY--XYY)«.25 
SNAP e 1 2 ) * ( 1 . " Y - X X + X X Y ) « . 5 

BHAPE(3)“-("'l . “XY+XX+YY-XXY+XYY)t«. 25 

SHAPE ( 4 ) = (1 . ♦-X ~YY -X YY ) . 5 

SIHAPE (.5) = I -1 . +XY+XX+YY+XXY+XYY)«. 25 

SHAPE<6)-'=( 1 . +Y-XX~XXY)».5 

SHAPE (7) "H, ”1 . -XY+XX+YY+XXY-XYYIw. 25 

SHAPE(8)=( 1 . -X-YY+XYY) W.5 . 

DERI V( 1 , 1 )■“ (Y+X2-XY2”YY)^is-25 
DERIV(1,2)=-X+XY ' 

DER I V ( 1 , 3 3 “ ( --Y+X2 - XY2+YY 3^.25 
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DERIV(1,4)!=(1. -YY)w.5 

DER I V (1 , 5 ) - ( Y+X2+X Y2+YY ) « . 75 

DERIVf. 1,6.^=-X"XY 

, DE: R I V (. 1 , 7 ) f. - Y+ X 2+X Y2 - Y Y ) » . 25 ' 

DERIV(l,8)=(-l.fYY)w.5 

I.)FRIV(2, 1) = <.X+Y2-XX-XY2)«.25 

DERIV(2,2.^=(-l.+-XX)w.5 

DERIV(2,3)=^R.-X+Y2-XX+XY2)s.25 

DERIV(2, 4.W--Y-XY 

DER I V ( 2 , 5 ) = (. X 4 Y2+ X X 4 X Y2 ) w . 25 

DERIV(2,6) = (1. -XX)w.5 

D E: f U V ( 2 , 7 > = C - X 4 Y 24 X X - X Y 2 ) « . 25 

DERIV(2,8) = "Y4-XY 

RETURN 

END 

C ■ 

BUBRUUT I ME D JACOB ( DER I V , DET JB , D J AC I , D JACK , I ELEM , 

! MELEM , MR 0 1 N , NNODP , LGAUS .) 

1 MPL I C I T REAL«8 ( A-H , 0 -7. ) 

DINENSION DERrV( 2 , 8 ),DJACI( 2 , 2 ),DJACK( 2 , 2 ) 

C0MMGN/ARC2/LN0DS(1 1 1, 8) , C00RD(386, 2) rNADEM(.386) , NDDFM(386) 
C SETUP TEMP MX TO ALLOW JACOBI IAN TO BE FORMED 
DO IDIME>1,2 
DO JDIME=1,2 
TEMPY=^0. 

DO INODP=l,NNODP 

KP 0 1 N- 1 ABS ( LNODS ( I ELEM , I NODP ) ) 

C KAG^'=(LGAUS-l)j<NNaDP4-INODP 

X X=DER I V ( I D I ME , I NODP ) wCOORD ( KPO I N , JD I ME ) 

TEMPY=TEMPY4-XX 

ENDDO 

DJACK ( IDIME, JDIME)=TEMPY 

ENDDO 

ENDDO 

DET JB^D JACK ( 1 , 1 )^^DJACK (2,2) -DJACK ( 1 , 2)sDJACK (2, 1 ) 

C CHECK FOR NEG OR ZERO DET 
If ( DET JB) 30, 30, 40 
30 CONTINUE 

TYRF.w, DET JB, •'DETJB' 

WR ITE ( w , 2000 ) I ELEM , LGAUS 

T YPEw, DJACK ( 1 , 1 ) , DJACK (2, 2) , DJACK ( 1 , 2) , DJACK (2, 1 ) 

STOP 

40 CONTINUE 

C INVERSE TEMP MX TO FORM JACOBIAN 
DJACK 1, 1)=DJACK(2,2)/DETJB 
DJACI (2, 2) =DJACK (1,1) /DETJB 
DJACI (1,2) “'’-DJACK! 1,2) /DETJB 
DJACI (2, 1 )="D JACK (2, 1) /DETJB 

2000 F0RMAT(//37H NON POSITIVE DETERMINENT FOR, ELEMENT, 214) 

RETURN 

END 

C 

SUBROUTINE SHAPE4(DERIV, SHAPE, XEOIV,YEDIV) 


IMPLICIT REAL?^8(A--H,0-Z) 

I.) IMENB ION DER IV (2 , 4 ) , BHAPE f 4 ) 

CwwLIN ELEMENT ANTICL NODE NUMBERING 
X^^-XEGIV 
Y-=V'EQIV 
XY^^X;«Y 

G1 lAPF ( 1) = U . -X -V+-X Y ) « . 25 
SHAPE (2) =f. 1. +X-V-XY)w.75 
SHAPE (. 3) f. 1 . +X+-Y+-X Y ) . 25 
3I1APE(4) -f. 1. -Xf-V~XY>w. 25 
DERIVC 1, 1) -<"1. +Y)«.25 
DERIU(1,2)=<.1.-Y)«.25 
DERIVCl, 3> = <'.1.+Y)«.25 
DERIU( 1, 4) =<-1. ~Y)w.25 
l)ERIU(.2, 1)=:<-1.+X)w.25 
DERIV(2, 2) =C-1. -X)w.25 
I)ERIU(2, 3)=C1. + X)w.25 
DERIV(2, 4> =U.-X)w.25 
RETURN 
END 

c 

SUE;« ROUTINE EKIUCOMIKBEL, KSURF,KBIDE., F'LUMX , EQRHS, VARB2, 
!MELEM,MPOIN, NNODP ,MTQTV, I ITER) 

IMPLICIT REALw8(A-H,0-Z) 

COMMON/ARC l/POSGP (3) , WEI6P (3) 

COMMON/ ARC2/LN0DB( 111,8), COORD (.386, 2) , NADFM (. 386 ), MODEM (386) 
COMMON /ARC 1 2/FR IFAC , P I , EN , COSA , S I NA , S I N2A 
COMMON/ARCia/RHO, MBOUEL 
CaMN0N/ARC30/ANGLN(386) 

CaMM0N/ARC14/IBEL(.60) , ISURFC60.I , ISIDE(60) 

DIMENSION EQRHStZOOO) , VARB 2 ( 2000 ) ,FLUMX(2a,28) 

D I MENS ION M C 3 ) , NODE ( 3 ) , DER I V ( 2 , 8 ) , SHAP 1(8), SHAPE ( 8 ) 
DIMENSION RH3U(8) , RHSV(8) , RHST(8) , CART(2, 8) , DERX (2, 8) 
DIMENSION C JACK ( 2 , 2 ) , C J AC I ( 2 , 2 ) , C05N ( 3 ) , B I NN ( 3 ) 

COMMON /ARC99 1 /TREE , PECLET , CONRAT , TFSTAR , HSTAR , 

!B5TAR,BIGY0, SPHTp,THKO,SIGMA 
CnMMON/AR(:887/SLETHG( 111,3) 

IDENTIFYING LOCAL NODE NOS 
IF(HSURF.LE. 3.)GaTO 88 
M(1 )-i 
M(2)^5 
M(3)-R 

NODE( l)"2sKSIDE-l 
NaDE<2)"'NODE(l,)+l 
N(JDE(.3)“MDDri(2.^ ’M 

I F ( MODE ( 3 ) . RT - B ) NODE ( 3 ) ^NODE ( 3 ) --8 
\ DO 145 ,S 

RHSUr. I) = 0. 0 
RHSV( I )-0. 0 
RHST( I)=---0, 0 
145 CONTINUE 

C GAUSS POINT EVAL STARTS 


/ 


LGAUS=0 

DO 150 IGAUS=1,3 

XEQIV=P0SGP(IGAUS3 

YEQIV=-1.0 

CALL SHAPES ( DERX, SHAP 1, XED^■;,YEQIv^^ 
TEMPR=0. 

DO :nodp=i,3 

KPG I N= I ASS ( LNODS K 3EL , I NGDP ; '■ 
ITrjTU=MADFM (KPQIN) 

I raTT=ITOTU+NODFM(KPOIN) -1 
EMPR=TEMPR-hVARB2 ITQtt • :^SHAF I •' INCDP • 
ENDDO 


C 

r 


CA!. 


PROP (TEMPR, SIGBTAR, CPSTAR, THKSTAR, 


DO KIMP=t,a 

KPOIN=IABSCLMaDS(KBEL, KINP ' 
rraTU=NADFh(KPOIN) 


I TGTT- 1 TOTUH-MODFM i . KPO T N ) - 1 
LGAUS=LGAU3+1 

GOTO C 1 0 , 20 , 30 . 40 ) , KS I DE 
CAL LENGTH OF ELcM AT EACH '3AUS3 PT 
10 CONTINUE 

XEQIV-'^POSGP^, [GAU3) 

VEaiV=-1.0 
ITENP^G 
SO TO CO 
20 CONTINUE 

XEaiV-=l . 0 

YEQIV^^A^aSGP '■ IGAUS) 

I TEMP =2 
GOTO SO 
30 CONTINUE 


XEQIV=POSGP (IGAUS) 
YEQIV=l.O 
ITEMP=L 


GOTO 50 
40 CONTINUE 


XEaiV=-l .0 
YEQIV==PnS6P( IGAUS) 


SIGN AY) 


j 7 Flvip =:2 
50 CONTINUE 

CALL SHAPE3CDERX,SHAPU XEQIV.YEQIV) 
CALL SHAPE4(DERIV, SHAPE, XEQIV,YEQrv) 


IELEM=KBEL 


MP0IN=:336 

MELEM«U1 


NN0DP=8 

CALL D JACOB (.DERX,DETJB,CJACI.C JACK, I ELEN, 

! HELEN , MPO I N , NNQDP , LGAUS ) 

TEMPY=DSORT ( CJACK ( ITENP , 1 ) ws2+CJACK UTENP , 2 )a#2) 


SLETH=TENPYwWei6P ( IGAUS) 
SLETHG(KBEL, IGAUS) =SLETH 
CAL Vise STRAIN RATE AT GAUSS PT FOR 


5 TH SURFACE ONLY 
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IF(KSURF.NE.5) GOTO 15 
DO 51 IDIME=1,2 
DO 51 INODP=l,NNODP 
CARTaDIME, IN0DP)=0.0 
DO 51 JDIME=1,2 

CART (IDIME, INODP)=CARTCIDIME, INODP .l+DERX C JDIME, TNODP > 
i^CJACKIDIME, JDIME) ' 

51 CONTINUE 

UVELY^O - 0 
VVELY=0. 0 
DUDX=0. 0 
DUDY=0.0 
DVDX=0.0 
DVDY=0.0 
DO 60 

KPOIN"=IABStLNODS(.KBELr I) ) 

ITaTU=NADFM(KPaiN> 

I TOT I TOTU+NODFM ( KPO I N ) -^2 
SHAP^=3HAP1( I) 

UVELY=UVELY+VARB2 ( ITOTU ) ^<^SHAP 
VVELY'=VVELY+VARB2 ( I TOTV ) -sSHAP 
CARXI”CART( 1, I) 

CARYI=CART(2r I) 

UVEL^=VARB2( ITOTU) 

VVEL~VARB2( ITQTV) 

DUDX=DUDX+CARXI«UVEL 
DUDY“DUDY<-CARYIi*UVEL 
DVD X --DVD X -("CAR X I « VVEL. 

DVDY=DVDY-+-CARY I wVVEL 
60 CONTINUE 

EPS^DSQRT •' 2 . Ow ( DUDX . 5w f. DUDY-hDVDX ) ?iw2-^DVDY;^^^2 ) ) 

I F ( EPS , LT .0.001) EPS==0 . 00 1 

VISC.= (9IGSTAR+B3TARwEPB^iw( 1 . 0/EN) ) / (1 .722;*EPS) 

V T = V VEL Y wCOS A -NJ VEL Y wS I N A 
C 

C EVAL GAUSS PT INTEGRALS FOR LOCAL RH5 VECTOR 
00 13 1=1,3 
NOD=NODE ( I ) 

SHAP=SHAPt(NOD) 

F ACTQR=^FR I F AC,»-V I SCwEPS 

RHSV ( NOD ) =RHSV ( MOD) -SHAP wFACTORwSLETH 

RHST(NOD ) =RHST (NOD)-i-SHAP’'sPECLET«FACTORiiVTwCQNRATT^SLETH/THKSTAR 
RHSU(N0D)=0. 

13 CONTINUE 

SKIPPING CONV BOUN CON FOR 5 TH SUR 
IF(KSURF. £Q.5)GOTO 150 
15 CONTINUE 

C IMPOSING SHEAR STRESS BC FOR 6 TH SUR 
IF(KSURF.EQ.6)THEN 
IP0INA=116 

IP0INB=122 I 

SL0PE=FACT0R/ (COORD ( IPOINB, 1 ) -COORD ( IPOINA, 1 ) ) 

C 

DO IJ=1,3 
NODD=NODE(IJ) 




SHAP =SHAP 1 ( NODD ) 

IPOIN=IABS(LNODB(KBEL,NQDD) ) 

XDIST=COORD( IP0IN8, 1) -COORD( IPOIN, 1 ) 

FACTOR=SLOPE«XDIST 

RHSV C NODD > =RHSV ( NODD ) -SHAP »FACTOR»SLETH 

RHST ( NODD ) =RHST CNDDD ) +SHAPsPECLETi^FACTORsVT:=fCQNRAT»SLETH/THKSTAR 
RHSU(N0DD)=0. 

ENDDO 

ENDIF 

C 

INCDR APPR CONV HT TR BC 
DO 75 ri = U3 
IN0DP=2»<KSIDE-1)+II 
IRQWU=7s(KSIDE-l)+M(I I) 

IFUNODP.GT.S) INDDP = INODP-a 
IF UROWU. 01.23) IR0WU=IR0WU-28 
KPOIN=IABSCLNODS(KBEL, INODP) ) 
rRaWT=rROWU+NODFM(KPOrN ) -1 
SHAP I =SHAP1( INODP) 

RHST ( INODP )=RHSTC INODP )+SHAPI^«HSTAR-o^TFSTAR-^iSLETH/THKSTAR 

DO 95 JC0N1=1,4 

JCOLP= ( JCDN 1 - 1 ) w7+2 

DO 95 JC0N2=1,2 

JNODP = ( JCON 1 - 1 ) w2+ JC0N2 

JCOLU= ( JCON 1 “ 1 ) »7+ JC0N2»4 -3 

JCGLT= ( JCON 1 - 1 ) w7+ JC0N2-ti3+ 1 

SHAP J=SHAP 1 (JNODP ) 

TERM“HSTARwSHAPI;o:SHAPJwSLETH/THKSTAR 
FLUMX ( IROWT, JCOLT) =FLUMX ( IROWT, JCOLT) +TERM 
95 CONTINUE 
75 CONTINUE 
ENDDO 

150 CONTINUE 
26 CONTINUE 

C INCOR NORMAL VEL BDC 
ANBLNC 100) =90.0 
DO 30 1 1 = 1 r 3 

C,«.^ EVAL DIR OF NORMAL AT EACH NODE 
GOTOCIOO, no, 120, 130) ,K3IDE 
100 CONTINUE 

XEQIV=FLOATCII) -2.0 
YEQIV=-1.0 
ITEMP=1 
RTEMP=+1.0 
SOTO 140 
no CONTINUE 

XEQIV=1.0 

YEaiV=FLOAT(II) -2.0 

ITEMP=2 

RTEMP=1.0 

GOTO 140 , 

120 CONTINUE 

XEQIV»-FL0AT(II)+2.0 ! 

YEaiV=1.0 

ITEMP=t 
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RTEMP=-1.0 
GOTO 140 
130 CONTINUE 

XEQIV=-1.0 
YEQ I V* -FLOAT ( 1 1 } +2 . 0 
ITEMP=2 
RTEMP=-1.0 
140 CONTINUE 

CALL SHAFE8CDERX,SHAP1,XEQIV,YEQIV) " 

EVAL slope at EACH NODE 
IN0DP=2^i(KSIDE-l)+II 
IRQWU=7»(KSIDE-1)+M(II) 

IFCIN0DP.GT.3) INODP=INODP -8 
IF(IR0WU.GT.28) IR0WU=IRQWU-23 
KPOIN=IABS(LNODStKBEL, INODP) ) 
IROWV=IROWU+NODFMCKPOIN) -2 
IR0WT=IR0WV+1 

IF(KSURF.EQ.5) ANGLN ( 100) =RAKE 
ANBLE= AN6LN ( KPQI N ) 

ANGLE=ANGLE»PI/180. 0 
C0SLX=RTEMP»DC0S C ANGLE ) 

COSN(II)=COSLX 

IF(DABS(COSN(II) ) .LT. l.OE-4) C0SNCII)=0.0 
COSL Y= -RTEMP DS I N ( ANGLE ) 

SINN(II)=COSLY 

IFCDABSCSINN(II)) .LT. l.OE-4) SINN(II)=0.0 

DO 90 JC0Nl=lr4 

JCOLP= t JCONl -1 ) «7+2 

DO 90 JC0N2=1,2 

JNODP= C JCON 1 - 1 ) »2+ JC0N2 

JCOLU= ( JCON 1 - 1 ) »7+ JC0N2»4 -3 

JC0LV=JC0LU+3-JC0N2 

JCOLT= JCOLU+4 - JC0N2 

SHAPPJ=SHAP1 CJNODP) 

CwwDET APPR MOM EON FOR NORMAL VEL BC 

IF(DABSCCQSNf,II) ) .GE.DABSCSINNdl) )i GOTO 55 
FLUMX ( I ROWV , JCOLU ) =COSN ( 1 1 ) wSHAPP J 
FLUM X C I ROW V , JCOLV ) =8 1 NN C 1 1 ) wSH APP J 
FLUMX ( I ROWV , JCOLT ) =0 . 0 
IF(JC0N2.EQ.2) GOTO 56 
FLUMX ( I ROWV, JC0LP)=0. 0 
56 RHSVC INODP) =0.0 

GOTO 90 
rnWTTIMUF 

FLUM X ( I ROWU , JCOLU ) =COSN C 1 1 ) wSH APP J 

FLUMX CIROWU,JCOLV)=SINN( II )»SHAPPJ 
FLUMX (I ROWU, JCOLT) =0.0 -- 

IF(JC0N2.Ea.2) GOTO 67 

FLUMX C I ROWU , JCDLP )=0 . 0 

67 RHSU( INODP) =0.0 
90 CONTINUE 
80 CONTINUE 

C«» ASS GAUSS PT INTEGRALS IN GLOBAL RHS VECTOR I 
DO 200 11=1,3 
NOD=NODE (ID 
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KPO I N'^ I ABS (LNODS (KBEL, NOD) ) 

' rTC:)H>NAl)FM(KPOIN) 

I TO TV™ I TO rU t-NODFM ( KPO I N ) -2 
ITOTT-"ITC:n U vNDDFM (KPOIN) -1 
EQRHSU TOTU)==EQRHS( ITOTU)+RHSU(NOD) 
EQf U<8 ( I T OT V ) -EQRHB ( I TOTV ) +RHSV ( NOD ) 
EORHS ( I TOTT ) =EQRHS ( I TOTT ) +RHST (NOD ) 
200 CONTINUE 
88 RETURN 
END 






C 

SUBROUTINE FORCE 


IMPLICIT realms ( A -H,0-Z) 

DIMENSION DERIV(2,8),DJACI(2,2),DJACK(2,2),SHAPE(8) 

D I MENS I ON CARTP (2,8), AREAW ( 9 ) , CART (2,10) 

D I MENS I ON C J ACK ( 2 , 2 ) , C J AC I ( 2 , 2 ) , COSN ( 3 ) , S I NN ( 3 ) 

DIMENSION M(3) ,N0DE(3) ,SHAP1(8) 

C0MM0N/ARC2/LN0DS( 1 11,8), COORD (386, 2) ,NADFM(386) ,N0DFM(386) 
C0MM0N/ARC13/RH0, MBOUEL 
COMMON/ARCaai /FW (111,3) 

COMMON/ ARC 1 0/CUTVEL , DCUT , W I D , RAKE 
COMMON / ARC 1 2 / FR I F AC , P I , EN , COS A , B I N A , B I N2 A 
C0MM0N/ARC3/VARB1 (2000) , VARB2(2000) 

COMMON/ARCl /POBBP (3) , WEIGP (3) 


COMMON/ AREA3/CARPQ(2, 72) , SHAP(3(72) , CARL6.(2, 36) , SHALG(36) 
C0MMaN/ARC991/TREF,PECLET,C0NRAT,TFBTAR,HBTAR, 

! BSTAR , S I GYO , SPHTO , THKO , S I OMA 


SET UP POSITIONS Z-. WTS FOR 3 POINT GUABB RULE 


POBGP ( 1 ) =0.7745966692 

P0SGP(2)=0.0 

P0SBP(3)=“P0BGP(1) 

WE I GP ( 1 ) =0. 5555555556 

WE I GP ( 2 ) =0 . 8888888889 

WEIGP(3)=WEIGP(1) 

KELEM=33 

TF0RCEL=0. 

FF0RCEL=0. 

C0NLEN=0 , 

XFGR=0. 

YF0R=0. 

XFC:)R1=0. 

YF0R1=0. 

DO 77 IE=1,8 

KELEM=KELEM+1 

L6AUB=0 

NGAUS=3 

TF0RCE=0. 

FF0RCE=0. 

F=0. 

FF=0. 

FNN=0. 
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non 


c CALCULATIONS START FOR EACH G.PT. 

DO 50 IGAUB=^1,NGAUB 
C SETTING UP GAUS POINT POSITIONS.' 

LGAUB=LGAUf3+l 
XEQIV=POSGP( IGAUS) 

ye:guv=-i.o 

C 

c CALCULATING PRESSURE. FOR G.PT. 

CALL SHAPE4(DERIV, SHAPE, XEQIV,YEQIV) 

TP=0. 

DO 40 INnDL=l,4 

IN0DP=aNaDL“l.)«2+l 

K PO I N~- 1 ABS ( LNODS ( KELEM , I NODP .) .) 

I0DFM=N0DFM(KP01N) 

IADFM:=NADFM(KPOIN) 

P--VARB1 (lADFM+l.^ 

TP •.•rip +.p „SH APE ( I NODL ) 

40 CONTINUE 
C 

CALL SHAPFeiDERIV, SHAPE, XEOI V, YEOI V.) 

lELEM^KELEM 

MF.LEM^lll 

MP0IN=386 

NN0DP==8 

CALL D JACOB (. DER I V , DET JB , C JAC I , C J ACK , I ELEN , 

! HELEN , MPO I N , NNODP , LGAUS 
C FINDING TEMPR. Z-. CALLING PROP VAR. 

TEMPR~0. 

' XBAUB=0. I 

YGAUS=0. 

DO 60 I NODP =1,8 

KPO I N= I ABS ( LNODS (HELEN , I NODP ) ) 

ITOTU=NADFN(KPOIN) 
rTaTT=ITOTU+NaDFN(KPOIN.l -1 
TENPR=TENPR+VARB2 ( ITOTT ) aSHAPE ( I NODP ) 

C FINDIND GLOBAL POS. OF G.PTS 

X GAUS= XG AUS+ ( SHAPE C I NODP ) ^^COORD ( KPO I N , 1 .5 .1 sDCUTs 1 000 . 

YG AU3=YQAUS+ ( SHAPE C INODP ) sCOORD ( KPO I N , 2 wDCUTw 1 000 . 

60 CONTINUE 

C CALLING PROP. VAR FOR TOTL G PT. TEMPR 

CALL PROP ( TENPR , S I GST AR , CPST AR , THKSTAR , S I GMAY ) 

C 

DO 51 IDIME=1,2 
DO 51 IN0DP=1,8 
I CART (IDINE, INODP) =0.0 

DO 51 JDIME=1,2 

C Af IT ( 1 1) I ME , I NODP ) =C ART ( I D I ME , I NODP ) +DER I V ( JD I ME , I NODP ) 
•wCJACKIDIME, JDINE) 

51 CONTINUE 

FINDING VIBC.,DUDX ETC., 

DUD X =0.0 
DUDY=0.0 
' DVDX=0.0 
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DVDY==0.0 
UVELY=^’0 . 

VVELY-O. 

DO 91 INODP=l,B 
KPOIN==IABS(L.NODS(KeLEM, INODP) ) 
ITaTU^=NADFM(KPOIN) 

I TOTV- 1 TOTU f NODFM t KPQI N) -2 
CAF^XI=CART(1, IMODP) 

CARYI=CART(2, INODPJ 
UVEL-VARB2(rTnTU) 

, VVEL=VARB2(IT0TV) 

UVEL.Y=UVELY+VARB2 ( ITOTU ) wSHAPE ( INODP ) 
VVEL Y“ VVEL Y+V/ARB2 f. I TOTV > wSH AP E ( I NODP ) 
DUDX=^ DUDX+CARX IwUVEL 
DUDY=DUDY+C ARY I wUVEL 
DVDX“DVDX+CARXIwVVEL 
DVDY=DVDY+CARY I wVVEL 
91 CONTINUE 


EPS=:T:)SQRT (2.0m (DUDX«w2+0. 5« (DUDY+DVDX ) mm2+PVDYws2) ) 
IF(EPS.LT. 0.0001) EPS=0.0001 

V I SO ( S I GST AR+BST ARmEPSmm ( 1 . 0/EN ) ) / (. 1 . 732mEPS ) 

VT-WELYmCOSA+UVELYmSINA 

FW (HELEN, IGAUS)=FRIFAC«VISCmEPSmVT 


C 

C NOW CAL. STRESSES FORCES USING 6.PT VAL.UES 
S I (3X X ■■•TP -1-2 . mV ISCmDUDX 
SIBYY=!"TP+2. mVISCmDVDY 
SI6XY^=VISCM(DUDY^fDVDX) 

F X - ( ■■■•?> I G X X mCOS A^ f S I GX YmS I N A ) 

FY=" ( -S I GX YmCOSA^<-S I GYY sS INA) 

IPDIN1 = IABS(LN0DS(KELEM, 1 )) 

I PO I N2= I ABS (LNODS (HELEN, 3 ) ) 

DE;L= DBQRT ( ( COORD ( I PO I N2 , 1 ) -COORD ( I PO I N 1 , 1 ) ) mm 2+ 
! (COORD( IP0IN2, 2) -COORD( IPOINt ,2) >sm2) 

DEL“I)EL/2. 

T FORCE-TFORCE t (FXmDELmWEI GP ( I BAUS) mWI D/DCUT ) 
FFORCE=FFORCE^»- ( FYmDELmWEI GP ( I GAUS) mWI D/DCUT ) 


C 

C CAL 

C CAL 


388 

50 


FRICTION Si NOF^NAL FORCES 
p~P4.FR I F'ACmV I SCmEPSmDELmWE I GP ( I GAUS ) mW I D/DCU T 
F ACTOR=FR I F ACm V I SCmEP B 
NORNAL FORCE FROM DERIVED REL 

STN0RI'i""(SIGXXwC0BAMM2-+'SIGYYMBINAMM2+SIGXYMSIN2A) 

FNN-FNN't-STNORNMDEI-MWE IGP ( I GAUS ) mWI D/DCUT 
ST SHE A= f. ( S I G Y Y -B I GX X ) /2 . mB I N2A-+'S IGXYm(1.-2.mSI NAmm 2 ) ) 
FF^FF^t-STSHEAwDELMWE I GP ( IGAU3.) sWI D/DCUT 
STNORMD-STNORMmSIGYO 


STSHE AD=STSHEAmS I GYO 

WR I TE ( 38 , M ) 'HELEN I GAUS X -C 

wRi reoa, m) #»### 


Y-C NORM. STRESS SHEAR STRESS-' 
############ ##########■' 


WRITE (38, 38B)KELEM, IGAUS, XGAUS, YGAUS, STNDRMD, FACTOR 
RMAT(/I3, 4X, 12, 2X, D8. 2, X, D8.2, D12. 4,2X,D12. 4.) 


FORMA 
CONTINUE 

XFOR~XFOR-+( 


-FmSINA-FNNmCOSA) 


I 
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YFOR=YFOR + ( -FwCOSA+FNNstSINA) 
XFOR 1 = X FOR 1 + (. -FFsS I N A -FNNsCOSA ) 
YFOR 1 = YFOR 1 + ( -FF»C0SA+FN1M»S I IMA ) 


C 


C 


C 


TFORCEL=TFORCEL+TFORCE 
FFORCEL=FFORCEL+FFORCE 
77 CONTINUE 

XFOR=XFOR^iB I feYOsDCUT^^DCUT 
YFOR=YFORwS IGYO«DCUT»DCUT 
XFOR 1 =XFOR 1 I GYO»DCUT«DCUT 
I YFORl=YFORl»SIQYO»DCUT»DCUT 

TFORCEL.=TFORCEL«S I GYOwDCUTkDCUT 
FFORCEL=FFORCEL»S I GYO»DCUT»DCUT 


199 


WRITE (35, 199) TFORCEL, FFORCEL, XFOR, YFOR, XFOR 1 , YFOR 1 


FORMAT (//8X, 'TANGENTIAL FORCE IS 


/8X, 

/QX, 

/8X, 

/8X, 


•'FEED FORCE 
FANG. FORCE 
■'FEED FORCE 
TANG. FORCE 


/8X,^FEED FORCE 


IS 

IS (FROM FS:N) 
IS (FROM FS.:N) 

IS (FROM FF?.:N) 
IS (FROM FFg.;N) 


■',D12.4,' N-', 
■■■,1)12.4,'' N-', 
■',D12.4,' N', 
■',D12.4,^ N-', 
■',D12.4,^' N', 
■',D12.4,-' N') 


RETURN 

END 




SUBROUTINE PROP (TEMPR, BIGSTAR, CPSTAR, THKSTAR,SIBMAY) 
IMPLICIT REALMa(A-H,0-Z) 

C0MM0N/ARC13/FIH0, MBOUEL, THKTOL 
COMMON/ARCIO/CUTVEL, DCUT, WID, RAKE, GAMMA 
C0MM0N/ARC12/FRIFAC, PI , EN, COSA, BINA, SIN2A 

COMMON/ ARC99 1 /TREF , PECLET , CONRAT , TFSTAR , H3TAR , 

! BSTAR, SIGYO, SPHTO, THKO, SIGMA 
TEMPR=TEMPRsTREF 
IF (TEMPR. GT. 373. ) GOTO 88 
SPHT=SPHTO 
SIGMAY=SIGYO 
THK=THKO 
GOTO 10 

88 IF ( (TEMPR. GT. 373. ) . AND. (TEMPR. LE. 1173. ) )THEN 
AM= -0.02226 
AC=50. 

ELSE 

IF (TEMPR. GT. 11 73) THEN 
AM=0. 00466666 
A 0^^29.63 
END IF 
END IF 

THK=AMsTEMPR^i-AC 

IF( (TEMPR. GT. 373. ) . AND. (TEMPR. LE. 473. ) )THEN 
0MX=0. 03185 
0CX==488. 119 

1*“ I 1*^" 

IF( (TEMPR. GT. 473. ) .AND. ( TEMPR. Le. 773. ) )THEN 
0MX=0.5859 


non 


0CX=229.239 
ELSE 

I F ' ( T EMP FT . 0T .773.). AND . ( TEMf^R . LE . 973 . ) ) THEN 
OMX--= 1.00435 
0CX=-~94. 222 
ELSE 

IF((TEMPR.GT.973. ).AND. (TEMPR. LE. 1 173. )) THEN 
OMX= -0.64865 
0CX=1514.146 
ELSE 

IF” ( (TEMPR. GT. 1173. ) .AND. (TEMPR. LE. 1373. ) )THEN 
0MX=0.0209 
0CX=728.76 
ELSE 

IF'( (TEMPR. GT. 1373. ) . AND. (TEMPR. LE. 1473. ) )THEN 
OMX-0. 4604 
0CX=^125.33 
END I F 
ENDIF 
ENDIF 
END IF 

ENDIF , 

ENDIF 

SPHT=OMX»TEMPR+OCX 
A^=2.5133E08 
B”191784.9 , 

CS== -366. 8211. 

A=9. 699 1731 EOS ' 

E= “2 137677. 

CS^=1521. ISO 

S I GMAY== ( A+B«TEMPR+CSwTEMPRss^2 ) 

10 SIGSTAR=SIGMAY/SIGYO 
CPEnAR=SPHT/SPHTO 
THKSTAR=THK/THKO 
RETURN 
END 

W W 

Ci<f WWiSW W WWfS S WWi'f'if SWSW WSS 

SUBROUTINE STRATES 
IMPLICIT REALw8(A-H,0-Z) 

D I MENS I ON DUD X ( 386 ) , DUDY ( 386 > , D VDX ( 386 ) , DVD Y ( 386 ) 
DIMENSION NCOUNT( 1000) , XPO) , YP(8) 

COMMON/ ARC992/STRATE (386) 

D I MENS ION DER I V ( 2 , 8 ) , D J AC 1(2, 2 ) , D J ACK (2,2), SFTAPE ( 8 ) 

D I MENS I ON C ARTP (2,8), ARE AW ( 9 ) , CART (2,1 0 ) 

D I MENS I ON C J ACK (2,2), C JAC 1(2,2), COSN ( 3 ) , 3 I NN ( 3 ) 

COMMON / ARC2 / LNODS (111,8), COORD ( 386 , 2 ) , N ADFM ( 386 ) , NODFM ( 386 ) 
COMMON/ ARC 1 2 / FR I FAC , P I , EN , COS A , S I N A , S I N2A 
C0MM0N/ARC3/VARB1 (2000) , VARB2(2000) 

COMMON / ARC 1 / P OSGP ( 3 ) , WE I GP ( 3 ) 

DATA XP/-1.0,0.,1.0,1.0,1.0,0.,-1.,-1./ 

DATA YP/-1.0, -1.0, -1.0,0. , 1. , 1. , 1. ,0./ 

NN0DP=8 
NELEM=1 1 1. 

C INITIALISING DUDX ETC. 
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DO INODP^l,NNODP 

DO IEL=1,NELEM 

I NOD- I ABS ( LNODS ( I EL , I NODP ) > 

DUI)X(IN0D)=^0.0 

DUDY(IN0D)=0.0 

DVDX(IN0D)=0,0 

DVDY( I NOD) =0.0 

NC0UNT(IN0D)=0 

ENDDQ 

ENDDO 


DO 200 IEL=1,NELEM 
DO 100 IN0DP=1,NN0DP 
I NOD= I ABS ( LNODS ( I EL, I NODP ) ) 

ITOTU=NADFM( INOD) 
rrOTV=ITOTUiMODFMC INOD) -2 
XEQIV=XP(INQDP) 

YEQIV=YP(INODP) 

CALL SHAPES ( DER I V , SHAPE , X EQ I V , YEQ I V ) 

IE.LEM=IEL 

MELEM=111 

MP0IN=386 

CALL D JACOB (DERI V, DETJB, CJACI , CJACK , lELEM, 
mELEM,riPOIN,NNODP) 

C 

DO 51 IDIME=1,2 
DO 51 INDP=1,8 
CART (I DIME, INDP)=0.0 
DO 51 JDIME=1,2 

CART ( I D I ME , I NDP ) =CART ( I D I ME , I NDP ) +DER I V ( JDI ME , I NDP ) 
IwCJACKIDIME, JDIME) 

51 CONTINUE 

CARXI=CART(1, INODP) 

CARYI=CART(2, I NODP) 

UVEL=VARB2(IT0TU) 

VVEL=VARB2CIT0TV) 

DUDX ( I NOD ) =DUDX ( I NOD) + CARX I ^t^UVEL 
DUD Y ( I NOD .) =DUDY ( I NOD ) +CARY I wUVEL 
DVD X ( I NOD ) =D VD X ( I NOD ) +CAR X I wV VEL 
DVDY( INOD)=DVDY( INOD) +CARYIwVVEL 
NCOl )MT ( 1 NOD ) ="NCOUNT ( I NOD ) + 1 
300 CONTINUE 
100 CONTINUE 
200 CONTINUE 
C 

DO IP0IN=1,386 

I DUD X ( I PO I N ) =DUDX (I PO I N ) /NCOUNI ( I PO IN) ^ 

DUD Y (I P 0 1 N ) = DUD Y ( I P 0 1 N ) / NCOUNT ( I P 0 1 N ) 

DVDX ( IP0IN)=DVDX ( IPOIN) /NCOUNT ( IPOIN) 

DVDYC IPOIN) =DVDY( IPOIN) /NCOUNT (IPOIN) 

C 

S I RATE ( IPOIN) =DSQRT (2. OwIDUDX ( IPOIN ) kw2 
! +0. 5=0! (DUDY( IPOIN) +DVDX (IPOIN) ) »m 2+DVDY( IPOIN) w»2) ) 
IF(STRATE( IPOIN) .LT. 0.0001) STRATE( IP0IN)=0. 0001 
WR ITE ( 90 , w ) I P 0 1 N , S TR A TE ( I P 0 1 N ) 
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ENDDO 

RETURN 

END 

SUBROUTINE CHTHK 
IMPLICIT REAL«8(A-H,0-Z) 

C0MM0N/ARC881 /FW( 1 11,3) 

COMMON/ARC12/FRIFAC, PI , EN, COSA, SINA, SIN2A 
COMMON/ARC13/RHO, MBOUEL 

COMMON / ARC 14/1 BEL f. 60 ) , I SURF ( 60 ) , I S I DE 1 60 ) 
COMMON/ARCIO/CUTVEL, DCUT, RAKE 

COMMON/ARC2/LNODS( 1 11,8), COORD (386,2) , NADFM (386) , N0DFM(.386) 
t COMMON/ ARC35/EPSN( 1 1 1 , 9) , XG ( 1 1 1 , 9) , YG ( 1 1 1 , 9) , VISCP (111,9) 

! ,DAREAG(111,9) ' 

C0MM0N/ARC887/SLETH6( 1 11,3) 

COMMON / ARC99 1 /TREF , PECLET , CONR AT , TFST AR , HST AR , 

! BSTAR, SI6Y0, SPHTO, THKO, SIGMA 
C 

TPLASTWK=0. 

DO KELEM=1, 111 
C 

PLASTWKG=0. 

DO KIGS=1,9 
C 

S I GFLO= 1 . 732«V I SCP ( KELEM , K I 6S ) sEPSN ( KELEM , K I GB ) 
PLASTWK=SIGFLawEPSN(KELEM,KIl3S) 
PLASTWKG=TPLASTWK+PLASTWK«DAREAG (KELEM, K IBS) 

ENDDO 

TPLASTWK=TPLASTWK+PLASTWKG 
' ENDDO 

TPLASTWKD=TPLASTWKi. ■ ' 

C' 

TFWK=0. 

1=33 

DO IJ=1,8 
I = I+-1 
FWKG=0. 

C 

DO IGAUS=1,3 

FWK=FW(.I, IGAUS)«SLETHG(I, IGAUS) 

FWKG=FWKG+FWK 

ENDDO 

TFWK=TFWK+FWKG 

ENDDO 

TFWKD=TFWK ’ . . ^ 

C 

T WKD=TPLASTWKD+TFWKD 

CHIPTK=(.COORD( 122, 1 ) -C00RD(386, 1 ) ) w;^2+ 

! ( COORD ( 386 , 2 ) -COORD (122,2)) ws2 

CHIPTK=DSQRT(CHIPTK)»DCUT«1000. 

^ WRITE (35, 223) CHIPTK , TPLASTWKD, TFWKD , TWKD ^ 

223 FORMAT!.//' CHIP THICKNESS = -,012. 4,' MM 

1 /•' TOTAL PLASTIC WORK = '',D12.4, ' ’ r 
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! /■' TOTAL FRICTION WORK =%D12.4, ' MM % 

! /■' TOTAL WORK = %D12.4, ' MM ' ) 

RETURN 
END 
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